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ABSTRACT 

A dual potential decomposition of the velocity field into a scalar and a vector 
potential function is extended to three dimensions and used in the finite-difference 
simulation of steady three-dimensional inviscid rotational flow and viscous flow. 

The finite-difference procedure has been used to simulate the flow through the 
80- by 120-Foot Wind Tunnel at NASA Ames Research Center. Rotational flow 
produced by the stagnation pressure drop across vanes and screens which are located 
at the entrance of the inlet is modeled using actuator disk theory. Results are 
presented for two different inlet vane and screen configurations. The numerical 
predictions are in good agreement with experimental data. 

The dual potential procedure has also been applied to calculate the viscous flow 
along two and three-dimensional troughs. Viscous effects are simulated by injecting 
vorticity which is computed from a boundary-layer algorithm. For attached flow 
over a three-dimensional trough, the present calculations are in good agreement with 
other numerical predictions. For separated flow, it is shown from a two-dimensional 
analysis that the boundary-layer approximation provides an accurate measure of the 
vorticity in regions close to the wall; whereas further away from the wall, caution 
has to be exercised in using the boundary-layer equations to supply vorticity to the 
dual potential formulation. 
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I. INTRODUCTION 
A. Overview of the Problem 

The numerical solution of the partial differential equations governing the conser- 
vation of mass, momentum, and energy is an important tool with which to analyze 
and study fluid motion. Recent advances in computer speed and memory together 
with improved methodology have resulted in computational fluid dynamics (CFD) 
playing an increasingly important role in aerodynamic research and development. 
The governing equations of fluid dynamics are the Navier-Stokes equations and sev- 
eral algorithms have been developed to solve these equations for simple two and 
three-dimensional configurations. However, even with recent developments in com- 
puter speed and memory, one of the pacing problems of today is the efficient and 
accurate prediction of flow 1 around complex three-dimensional geometries. 

Several approaches are commonly employed in order to reduce the overall com- 
putational time and cost of obtaining numerical solutions of steady-state flow prob- 
lems. One strategy is to solve a single set of simplified equations to approximate the 
entire flowfield. The most commonly used simplification for high Reynolds number 
flows is to solve the full potential equation, thereby assuming the flow to be inviscid 
and irrotational. However, even at high Reynolds number, many of the flows in 
nature possess substantial regions with rotation or vorticity. If one were to solve 
these problems with the potential equation, then the tedious task of setting up 
vortex-sheet discontinuities in the field and adjusting them to the surrounding flow 
has to be encountered. The alternative, which is the next step in the refinement of 
the approximation, is to solve the Euler equations which permit rotational flow ev- 
erywhere. Finally, if one is interested in the viscous effects at solid boundaries, then 
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it is necessary to solve the complete Navier-Stokes equations or an appropriately 
reduced form of these equations which contains viscous terms. 

Another strategy which is often employed is to divide the flow-field into separate 
zones and treat each zone with a simplified set of equations with matching performed 
in order to couple the separate solutions. The equation sets, and hence the solution 
algorithms, used in the different regions may be drastically different as. for example, 
in the commonly used viscous-inviscid interaction methods in which a boundary- 
layer algorithm is coupled with an inviscid flow solver. 

After having selected a simplified or complete equation set, one is faced with sev- 
eral alternative procedures for solving them. For two-dimensional flows, it is as com- 
mon lo find solution algorithms written to solve for the primitive variables: velocity 
and pressure, as it is easy to find algorithms to solve for the derived variables: stream 
function and vorticity. In two dimensions, the advantage of the vorticity/stream- 
function approach is that cont inuity is automatically satisfied at the expense of solv- 
ing a higher-order differential equation. Also, the same vorticity /stream-function 
approach results in fewer variables due to the elimination of the pressure from the 
governing equations. For the more general three-dimensional case, this approach in 
fact leads to more variables than the primitive variable case. However, the method 
may have advantages in certain applications. 

The focus of this study is to develop a procedure for solving the equations of 
motion with a vorticity /stream-function type of method, and a review of previous 
work done in this area is given in the next section. 
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B. Literature Review 

L Vector potential methods 

a. Incompr essible flow Several efforts have been made to extend the ad- 
vantages offered by the vorticity/stream-function approach in two dimensions to 
three-dimensional flow problems by expressing the velocity as the curl of a vec- 
tor potential so that continuity is automatically satisfied (for incompressible flow) 
since the divergence of the curl of a vector is identically zero. A review of these 
efforts is provided by Richardson and Cornish (1977) and by Morino (1985). Be- 
cause early efforts were, in part, hampered by lack of computational resources, not 
much work was reported for three-dimensional problems until Hirasaki and Hel- 
iums (1967) formulated the equations and the boundary conditions by expressing 
the velocity as the curl of a vector potential. Aziz and Heliums (1967) applied 
this formulation to study three-dimensional natural convection in enclosures. The 
formulation of Hirasaki and Heliums required the solution of a second-order partial 
differential equation in order to obtain the boundary conditions on the vector po- 
tential. Subsequently, they realized (Hirasaki and Heliums. 1970) that the addition 
of a scalar potential function to the velocity decomposition simplified the bound- 
ary conditions on the vector potential at the expense of introducing an additional 
variable and its corresponding equation. This velocity decomposition, in a slightly 
different form, has been termed the “dual potential” decomposition by Chaderjian 
and Steger (1983) and is used in this study to refer to the decomposition of the 
velocity into a scalar and a vector potential. The vector potential has been found 
to be particularly useful in the solution of three-dimensional natural convection in 
enclosures and several results have been reported ( cf . Mallinson and de Vahl Davis, 
1973; Ozoe et ah, 1976; 1977; 1982). For these problems, it is not necessary to 
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include the scalar potential and hence computation time and storage requirements 
are reduced. Aziz and Heliums (1967) also showed that the vector potential for- 
mulation is computationally less expensive than a comparable primitive variable 
formulation for the problems they considered. 

Aregbesola and Burley (1977) applied the dual potential formulation of Hi- 
rasaki and Heliums (1970) to study incompressible flow through ducts with inter- 
nal recesses. Richardson and Cornish (1977) presented a rigorous analysis for the 
existence and uniqueness of the vector potential for general three-dimensional in- 
compressible flow problems. In order to solve the full Navier-Stokes equations, the 
dual potential formulat ion requires the solut ion and storage of seven variables (four 
potential functions and three vorticity components) as opposed to four (three ve- 
locity components and pressure) in the primitive variable formulation. Hence, there 
have been several attempts to reduce the number of variables. Wong and Reizes 
(1984) replaced the scalar potential function with a specified inlet velocity vector, 
thereby reducing storage and computation time. However, their technique was re- 
stricted to flow through constant area ducts with incoming irrotational flow. More 
recently, Yang and Camarero (1986) reported solutions for three-dimensional duct 
flow problems using a dual potential type formulation. They used the Hirasaki and 
Heliums (1970) formulation for the boundary conditions on the vector potential. In 
particular, their formulation allowed for incoming rotational flow. 

In recent years, with the availability of greater computer resources, some of the 
advantages of using the vector potential formulation has revived interest in this ap- 
proach. One area of ambiguity in the boundary conditions for the vector potential 
has been the case of multiply connected regions such as annular passages. Richard- 
son and Cornish (1977) derived the boundary conditions on the vector potential for 
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general three-dimensional geometries including multiply connected regions. Wong 
and Reizes (1986) derived slightly different boundary conditions for such regions 
and presented solutions for flow through cylindrical annular passages. 

Related methods have been presented by Davis et al. (1986) who used the 
decomposition of the velocity field into two stream functions and applied it to 
compute the incompressible viscous flow over three-dimensional troughs and bumps. 
An excellent review of the techniques used for solving the Navier-Stokes equations 
in two dimensions with the vorticity/stream-function method is given by Roache 
(1972). 

b. Co mpressible flow All of the studies mentioned up to now have dealt with 
incompressible steady flow problems. In inviscid compressible flow, the potential 
equation has been the workhorse of the aerospace industry, being widely used to 
compute the flow around complex configurations for ?. wide range of flow conditions. 
In the presence of shock waves, the flow is no longer irrotational, and entropy and 
vorticity corrections are needed to give better estimates of flow parameters. Em- 
mons (1948) was the first to use the stream function to account for rotational effects 
and entropy changes following shocks in two-dimensional transonic flow. Hafez and 
Lovell (1981, 1983) also used the stream function to develop procedures to correct 
the potential solution. They studied a variety of methods to include rotational ef- 
fects in the potential solution. Chaderjian and Steger (1983, 1985), working along 
similar lines, developed the dual potential formulation to simulate the steady tran- 
sonic inviscid rotational flow around airfoils. Atkins and Hassan (1983) presented a 
different stream function formulation in which density and stream function were the 
dependent variables and the governing equations were solved in strong conservative 
law form. Recently, Hafez et al. (1987) developed a finite-element formulation in 
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which they used a dual potential type velocity decomposition scheme and applied 
it to calculate the transonic viscous flow through symmetric nozzles. 

In three dimensions, Chaviaropoulos et al. (1986) developed a method for 
solving compressible inviscid rotational flow in the subsonic regime and used it 
to simulate flow through a rectangular elbow. Their formulation is similar to the 
one presented in this work for inviscid flow. Other alternate formulations using 
the Clebsch velocity decomposition have been presented by Grossman (1983) to 
compute inviscid supersonic conical flows, and by Chang and Adamczyk (1983) to 
compute inviscid subsonic rotational flow in turning channels. Sherif and Hafez 
(1983) presented a scheme for solving irrotational transonic flow using two stream 
functions and discussed possible extensions to rotational flow. Morino (1985) de- 
rived a set of equations for unsteady, compressible viscous flow. No implementat ion 
of the formulation has yet been reported. 

2. Vi scous flo w 

A number of methods have been used to simulate three-dimensional viscous 
flows. Some of the recent developments and trends in viscous flow simulation have 
been reviewed by Steger and Van Dalsem (1985) and by Shang (1985). Viscous 
flow can be simulated either by solving the complete Navier-Stokes equations (or a 
reduced form of these equations) or, alternatively, by coupling the boundary-layer 
equations to an efficient inviscid flow model (e.g., Lock, 1981; Melnik, 1981; Mc- 
Donald and Briley, 1983). Various approximations to the Navier-Stokes equations 
such as the “parabolized’' Navier-Stokes equations (PNS), the “partially parabo- 
lized” Navier-Stokes equations (PPNS), and the “thin-layer” Navier-Stokes equa- 
tions (TLNS), are discussed by Anderson et al. (1984). 
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Efforts at solving the Navier-Stokes equations have been aimed at improving the 
robustness and speed of the various algorithms that are currently popular (Beam 
and Warming, 1976; Briley and McDonald, 1977; MacCormack, 1982). All of these 
algorithms solve the governing equations in primitive variables. Other algorithms 
have been developed for solving the incompressible Navier-Stokes equations ( cf . 
Orszag and Israeli, 1974). A survey of the present state of the art in the numerical 
solution of the Navier-Stokes equations is given by Holst (1987), while details of 
the various methods are given by Anderson et al. (1984). It is sufficient to mention 
here that the increased availability of supercomputers has made it possible to per- 
form three-dimensional calculations of complete aircraft configurations (Shang and 
Scherr, 1986; Flores et ah, 1987). However, due to time and storage limitations, it is 
st ill necessary to find ways to improve the convergence rate of the Navier-Stokes al- 
gorithms. One approach has been the “fortified" Navier-Stokes concept introduced 
by Van Ilalsem and Steger (1985). In this approach, the Navier-Stokes equations 
are solved in the entire flow field. In addition, the boundary-layer equations are 
solved near the wall on a fine grid to resolve the viscous gradients near the wall. 
Th is boundary-layer solution is then superimposed as a source function by inter- 
polation to the main grid. This approach was found to improve the convergence 
rate of the Navier-Stokes algorithm by a factor of 20 for the cases considered (Van 
Dalsem and Steger, 1986b). Along similar lines, Goble and Fung (1987) developed 
a truncation error scheme in which a local solution was obtained on a fine grid and 
used to form an approximation to the truncation error. This approximation was 
then used as a forcing function in the global solver. 

On the other hand, the search for alternative faster, but possibly not as general, 
methods has resulted in various zonal schemes. Here, efforts have been concentrated 
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on making the schemes more general by improving the approximations made in the 
various zones. One such approach to improve the accuracy of the inviscid flow 
approximation was proposed by Steger and Van Dalsem (1985) wherein viscous 
effects were accounted for by the injection of vorticity into an inviscid rotational 
solution. The vorticity was obtained from a suitable boundary-layer solution which 
interacted with the inviscid rotational flow through the pressure field. This ap- 
proach was implemented in two dimensions to calculate the separated flow through 
a diffuser. A similar idea, referred to as the “pseudo Navier-Stokes approach”, was 
proposed at about the same time by Whitfield (1985) with the difference that Steger 
and Van Dalsem used the vorticity /stream-function method to calculate the invis- 
cid flow while Whitfield maintained the primitive variable formulation for both the 
boundary-layer and the outer flow equations. Halim and Hafez (1984) developed 
a slightly different scheme in which they used the stream function to calculate the 
inviscid How and the “partially parabolized” Navier-Stokes equations to supply the 
vorticity instead of the boundary-layer equations. 

One very useful approach for resolving viscous effects adjacent to solid surfaces 
and wakes is the viscous-inviscid interaction method. In this method, the viscous 
shear layer over the wall or the wake, and the inviscid flow external to the shear layer 
are solved separately, together with a matching process which allows for the interac- 
tion of the two flow’s. In the most common viscous-inviscid interaction schemes, the 
viscous layer is resolved by solving the boundary-layer equations, and the inviscid 
flow by solving the full-potential equation. The two solutions are matched by iter- 
ating for the displacement thickness. Some like Whitfield et al. (1981) have used 
the Euler equations to obtain the inviscid flow solution. In two dimensions, vari- 
ous methods have been proposed to improve the speed and accuracy of the viscous 
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and the inviscid sets of equations ( cf . Lock, 1981; Melnik, 1981; McDonald and 
Briley, 1983). Most of these have concentrated on obtaining the quasi-simultaneous 
or simultaneous solution of the two sets of equations (Lee and Pletcher, 1986), 
to improve the convergence rate. Initially, viscous-inviscid interaction procedures 
used the direct method in which the pressure was calculated from the inviscid solu- 
tion and specified as a forcing funct ion for the boundary-layer equations. However, 
with the pressure specified, the boundary-layer equations become singular at sep- 
aration (Goldstein, 1948; Brown and Stewartson, 1969) and the equations must 
be solved in the inverse mode (Catherall and Mangier, 1966; Klineberg and Steger 
1974; Williams, 1975; Carter and Wornum, 1975; Kwon and Pletcher, 1979). In 
the inverse mode, the direct forcing function (pressure) is replaced by specifying 
an inverse forcing function such as displacement thickness or skin friction. In three 
dimensions, it is sufficient to only solve the streamwise momentum equation in the 
inverse mode to avoid the separation-point singularity, while the cross-flow momen- 
tum equation can be solved in the direct mode (Edwards and Carter, 1985; Van 
Dalsem and Steger, 1986a), though results have been obtained with both momen- 
tum equations solved in the inverse mode, that is, both components of skin friction 
or displacement thickness specified (Cousteix and Houdeville, 1981; Wigton and 
Yoshihara, 1983; Delery and Formery, 1983; Radwan and Lekoudis, 1984; Edwards 
and Carter, 1985). 

There have been a number of schemes proposed for discretizing the steady three- 
dimensional boundary-layer equations (Der and Raetz, 1962; Dwyer, 1968; Krause 
et al., 1976; Kitchens et al., 1975). Most of these studies address the problems asso- 
ciated with the hyperbolic character of the equations in the x-z plane. In addition, 
for steady flows, another problem is the requirement of initial data along a y-z plane. 
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An excellent review of the contributions in this area can be found in the book by 
Anderson et al. (1984). Most of these problems are circumvented by Van Dalsem 
and Steger (1985) by using a time-dependent approach for obtaining solutions for 
steady flows. Even though only steady-state solutions are computed, with the use 
of the unsteady equations, relaxation algorithms can be devised which avoid some 
of the limitations of standard space-marching schemes. In addition to being able 
to compute unsteady flows, solving the unsteady equations may provide a simpler 
and alternative method for computing steady, separating and reattaching flows by 
enabling the use of flow-dependent difference operators. Additionally, with the use 
of a relaxation approach, simpler boundary conditions can be specified. This is 
especially true for boundary-layers in three dimensions, where, in a space marching 
scheme, all variables must be specified at inflow and at least one side boundary; 
whereas in a relaxation scheme, simpler zero-gradient boundary conditions can be 
used at the side boundary and an additional set of reduced equations for a plane of 
symmetry need not be solved. 

C. Scope of the Present Study 

This study is a combination of several of the solution strategies outlined above. 
The goal of providing a more complete solution to the flow-field is achieved by using 
a zonal method wherein a simplified equation set - namely the potential equation - 
is employed in regions of irrotational flow which is coupled with the Euler equations 
in regions of rotational flow. In the present method, a simplified set of equations is 
derived which have several useful features: 

1) For the simplest kind of flow, without any modifications to the numerical 
algorithm, the method reduces to solving the full-potential equation. Hence, it is 
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possible to take advantage of the many fast algorithms that have been developed 
for this widely used equation. 

2) For inviscid rotational flow, it is necessary to solve additional linear Poisson 
equations to account for rotational effects. Again, fast solution schemes can be 
used. 

3) Finally, with the present approach, it is easy to extend the method to more 
general flow problems by solving additional equations for transport of vorticity or 
alternatively couple the inviscid set in a novel manner with a suitable boundary-layer 
algorithm, thereby retaining the advantages of the individual solution algorithms 
for the simplified sets of equations. 

In this study, the dual potential formulation of Chaderjian and Steger (1985) is 
extended to three dimensions and used to simulate the steady, inviscid rotational 
flow through indraft wind tunnels. Screens and vanes which are located at the 
inlet entrance are modeled using actuator disk theory. The numerical results are 
compared with experimental data. A new interaction scheme is developed for cou- 
pling the inviscid dual potential formulation with the three-dimensional boundary 
layer procedure developed by Van Dalsem and Steger (1985). This interaction pro- 
cedure is applied to compute the steady, attached flow over a three-dimensional 
trough. The procedure is also used to compute the steady, separated flow over a 
two-dimensional trough using a prescribed vorticity boundary condition. 

The conservation laws for steady, three-dimensional flow are presented in Chap- 
ter II. Also in Chapter II, the dual potential velocity decomposition scheme is pre- 
sented together with the numerical algorithms for solving the various equations. In 
Chapter III, the solution procedure for inviscid rotational flows is presented, along 
with results for flow through an indraft wind tunnel. This chapter also includes a 
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description of the actuator disk theory used to model the screens and vanes which 
are located at the entrance of the wind tunnel. Chapter IV presents the solution 
method for viscous flows along with results for flow over various trough configura- 
tions. This is followed by the conclusions in Chapter V. Appendix A includes an 
analysis of the geometric errors created by using nonuniform grids in transformed 
coordinates, while Appendix B gives the numerical algorithms for solving the in- 
compressible two-dimensional dual potential and boundary-layer equations. 



II. DUAL POTENTIAL FORMULATION 

The objective of the portion of the study described in this chapter is to de- 
velop a formulation for simulating steady, three-dimensional rotational flow using 
the dual potential velocity decomposition. With the dual potential decomposition, 
the dependent variables are derived variables, vorticity and entropy. In order to 
derive the necessary equations for both inviscid and viscous rotational flow, the 
governing equations of mass, momentum and energy will first be presented in the 
more conventional primitive variables where velocity components and pressure are 
the dependent variables. Then the dual potential form of the equations will be 
derived. Finally, the numerical algorithms for the various equations will be given. 

A. Conservation Laws 

1 . G ov e r n i ng equations in primit i ve vari ables 

The governing conservation equations of mass, momentum and energy for the 
steady, three-dimensional flow of a perfect gas are derived in many textbooks in 
fluid mechanics. Using the form given by Anderson et al. (1984), the equations are 
presented below in Cartesian coordinates. 

Continuity: 

Vpq = 0 (2.1) 

Momentum: 

p{q ’ V)g = -Vp + V - f (2.2) 


Energy: 


p(qV)H^V-(rq-f) 


(2.3) 
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Equation of State: 


p — pRT 


(2.4) 


In these equations, p is the fluid density, q — (u,v, w) T is the velocity vector, p is the 
pressure. T is temperature, H is the total enthalpy, R is the perfect gas constant, 
V is the gradient operator, and f is the viscous stress tensor given by 


‘ij 


/' 




2 du k 
3 t] dx k 


(2.5) 


where is the Kronecker delta (6 t j = 1 if i = j and S t j = 0 if t ^ j), and p is the 
viscosity. Using Fourier’s law for heat transfer by conduction, the heat flux vector, 
/ can be expressed as 

/ - - kXT (2.6) 


where k is thermal conductivity. The total enthalpy H is defined in terms of the 
enthalpy h as 

H=h+ ] 2 q-q (2.7) 

Also, for a perfect gas, the following relationships hold: 

h = c P T] c P -c v =R (2.8) 

c X) 

where c p is the specific heat at constant pressure, c v is the specific heat at constant 
volume and q is the ratio of specific heats. 


2_. Nondimensionalization 

Variables in the governing equations are first nondimensionalized by their free 
stream values. The nondimensionalization is shown below with the nondimension- 
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alized variables indicated by the carets (“) . 


X — 

X 

L 

y = 

y 

L 

Z — 

2 

L 


u 

v = 

V 

w — 

w 

u = 

u oc 

Uoc 

u oc 


P- Poo 


J>_ 


H 

p = 

P U lo 

P = 

P oc 

H = 

V > 
u oc 


In these equations, L is a reference length. In all subsequent 


M 


M = 

Moo 


(2.9) 


material, unless other- 


wise noted, variables will be assumed to be nondimensionaiized although the carets 


will be omitted for ease of reading. 


B. Velocity Decomposition Scheme 


It is well known that any velocity field can be expressed as the gradient of a 
scalar potential and the curl of a vector potential (c/. Aris, 1962; Panton, 1984): 

t? 


V<£ + V >■ B\ B 


( 2 . 10 ) 


so that 


bar T Ipy Xz 

q = ( v 1 « I <i>y + d z *- rp x 

w / \<t>z + Xx - tiy 

With this decomposition, the divergence of q is given by 

V • <f = 


( 2 . 11 ) 


( 2 . 12 ) 


where V 2 is the Laplacian operator V • V, and the curl of the velocity (or vorticity) 
is given by 

( Wy — V z \ ( ~(l?yy + $ zz) + ( Xxy T V’xz) \ 

u z — w x I = I ~{Xxx + Xzz ) + {$yx + 0j/z) I (2-13) 

Vx~ U y / \ — (ipxx + 4>yy) + {#xz + Xzy) ) 

A rigorous mathematical analysis for the existence and uniqueness of the vector 

potential functions was presented by Hirasaki and Heliums (1967), and by Richard- 
son and Cornish (1977), and proofs will not be given here. The three-dimensional 
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velocity field is represented by four potential functions, so there is one degree of 
redundancy in this decomposition. Hence, a constraint needs to be imposed on 
these functions in order to uniquely specify the velocity field. Following Chaderjian 
(1984), and Chaderjian and Steger (1985), the scalar potential is retained in order 
to obtain an equation for 4> that is similar to the transonic full potential equation 
and a constraint is therefore imposed on the vector potential functions. A standard 
vector potential constraint that has often been employed (cf. Panton, 1984) is the 
condition that the vector potential be divergence free, that is, 

V • B = d x + Xy + 4’z - 0 (2.14) 

Th is consistency expression relates the vector potential functions to one another 
and removes the redundancy. It is used both to simplify the vorticity relations and 
to obtain an extra boundary condition (see also Hirasaki and Heliums, 1970; Wong 
and Reizes, 1984: Chaviaropoulos et ah, 1986). 

Subject to the consistency condition of equation (2.14), the vorticity assumes 
the compact symmetric form: 

V 2 tf = - W] 

V 2 X = - W 2 (2.15) 

v 2 v> = - uj :i 

or simply 

V 2 B - -w (2.16) 

An additional consistency equation relating the vorticities is obtained from the 
vector identity 


V*u5 = V- Vx<f=0 


(2.17) 
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or 

( w l)x + (w-2)y + (^ 3)2 — 0 (2-18) 

which is also the Laplacian operating on equation (2.14). 

L Bo undary co ndit ions 

On any boundary surface two linear combinations of the vector potential func- 
tions can be kept constant, and the vector potential consistency relation, equation 
(2.14), can be used to solve for a third linear combination. Tangency is then imposed 
on a body boundary surface from 

n -V<i> = -n -V x B (2.19) 

In Cartesian coordinates, two components of B can be chosen to be constant on an 
T.y, or 2 boundary plane so that n • V x B - 0 and V • B — 0 supplies a Neumann 
boundary condition for the third component. For example, on a 2 = constant plane, 
and \ are taken as constants. Then and \ y — 0 and consistency gives ip z = 0, 
from which 0 can be determined on the 2 — constant surface. Likewise, \ x and 
fly ~ 0 and tangency, w = 0 , is satisfied by setting 4> z = 0 . 


2. The continuity equation 


For the velocity decomposition given by equation (2.10), the continuity equation 
(2.1) takes the form 

( P<t>x)x + ( P<t>y)y + {p<t>z]z = 


(2.20) 

— ~~ Xz)px “I" {&z ~ V 7 x)Py + (Xz — fly)Pz\ 

It may be noted that this differs from the usual form of the transonic full-potential 

equation by the presence of the nonzero right-hand-side terms describing the ro- 
tational component of the velocity. For constant density, equation (2.20) reduces 
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simply to 

V 2 <£ - 0 (2.21) 

It is useful to note that for a given vorticitv distribution, the dual potential set of 
equations given by equations (2.16) and (2.20), subject to ihe appropriate boundary 
conditions, can be solved to obtain the complete velocity field. In particular, for 
the vector potential functions, these equations are simple linear Poisson equations, 
which are weakly coupled and hence can be individually solved efficiently using 
existing well developed methods. The dual potential set of equations outlined above 
are obtained from the continuity equation and from mathematical simplifications 
developed from the velocity decomposition. For irrotational inviscid flow, vorticity 
is zero and only the full potential equation needs to be solved, as the momentum 
equations are, in effect, satisfied. For rotational flow, however, the momentum 
equations which prescribe the vorticity need to be considered in order to obtain a 
complete description of the flowfield. These equations are developed in the next 
section for both inviscid and viscous flow. 

C. The Momentum Equations 
1_. Momentum equations for inviscid flow 

Inviscid flow is governed by the Euler equations which are obtained from the 
governing equations (2.1) - (2.8) by neglecting viscous and heat conduction terms. 
After omitting the viscous terms, the momentum equation (2.2) can be rewritten 
as 


p{q- V)q = - Vp 


( 2 . 22 ) 
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Making use of the identity 

(2.23) 

equation (2.22) can be rewritten as 

v(£)-f*a=-^ (2.24) 

Crocco’s equation is obtained from the momentum equation by replacing the 
pressure gradient with an entropy gradient by using Gibb’s equation which is ob- 
tained from the First and Second Laws of Thermodynamics and can be written 
as 

TVs = Vh - — (2.25) 

P 

Crocco’s equation for a steady, adiabatic flow ( h -+ \q 2 = constant ) can then be 
w'ritten as 

<f x u; = -TVs (2.26) 

An equation for the convection of entropy is formed by taking the dot product 
of the velocity vector with Crocco’ equation: 


q ■ Vs = 0 


(2.27) 


Note that equation (2.27) is formed from a linear combination of the Crocco equa- 
tions, and thus replaces one of the three equations given by equation (2.26). 


Also for a steady inviscid adiabatic flow of a perfect gas, Crocco’s relations, to- 
gether with the perfect gets relation, can be used to derive the compressible Bernoulli 
equation (c/. Anderson et al., 1984): 


Pr 


. 1 ( 1 ) 

2\h r ) 


1 

1 ~ 1 


( S -Sr)/i? 


(2.28) 



where the subscript r refers to a reference state. 


These equations are modified for incompressible flow. For incompressible flow, 
density is taken to be constant and Bernoulli’s equation is no longer used. In 
addition, the total pressure, Po = P + \p<] 2 , is used in place of entropy, and Lamb’s 
equation replaces Crocco’s equation. Making use of equations (2.24) - (2.26), the 
momentum equation for inviscid incompressible flow is written as Lamb’s equation: 


<f x u> = 


Yi o 

p 


(2.29) 


and the entropy convection equation is replaced with a total pressure convection 
equation: 

q- Vp 0 - 0 (2.30) 


Note again that the total pressure convection equation was formed by taking the 
dot product of the velocity vector with Lamb’s equation and therefore replaces one 
of the three equations given by Lamb’s equation. 

Boundary conditions for the vorticity and entropy/total pressure are dependent 
on the geometry and flow conditions and will be given in the Chapter III, where 
results will be presented for a particular application of the procedure developed 
here. 


2. Momentum equations for viscous flow 

For viscous flow, the more complete momentum equations need to be solved. 
Assuming constant viscosity and density, equation (2.2) can be rewritten as 
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where Re is the reference Reynolds number defined as Re = PocUocL/ Poo- Taking 
the curl of equation (2.31) and simplifying gives 

p{q ■ V)w = (w • V)g + -^-V 2 w (2.32) 

The momentum equation in this form, together with the dual potential equa- 
tions can be solved to obtain the complete flowfield. However, in most cases, it is 
known that viscous effects are confined to a thin layer adjacent to the wall and most 
of the outer flow is inviscid. This is the basis for boundary-layer theory. The effect 
of the boundary-layer on the outer inviscid flow is usually accounted for through the 
use of a displacement effect (Lighthill, 1958). With the present formulation, it is 
possible to develop an alternate representation for the influence of viscous effects on 
the inviscid flow - that is, by obtaining the vorticity in the viscous region from the 
boundary-layer equations instead of solving the vorticity transport equations given 
by equation (2.32). That is, it may be possible to take advantage of the speed and 
efficiency of boundary-layer algorithms to calculate vorticity. With the vorticity 
specified, the dual potential set of equations could be solved to obtain the complete 
flow field. This concept is not new and has been proposed and implemented in 
two dimensions by Steger and Van Dalsem (1985). Various schemes for using the 
boundary-layer equations to supply the vorticity to the dual potential equations 
will be described in Chapter IV. 

D. Generalized Curvilinear Coordinates 

In order to simplify the treatment of arbitrary body boundaries in the numerical 
simulation, body-fitted curvilinear coordinates are employed, and the flow domain 
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is mapped to a uniformly spaced rectangular coordinate region. The general coor- 
dinate transformation is defined by 

£ = £{x,y,z), ri = ri(x,y,z), . f = f(x,y, z) (2.33) 

In terms of the independent variables the Cartesian velocity components can 

be evaluated from the potential functions using the chain rule of partial differenti- 
ation on equation (2.11) as 

u = (far 0c + T] X 4) V + + (ft + r/yt/'y + CyVt;) 

- (CaXf + ^Xy + fcX?) 

v = (f„«/>tf + ri^ + f v <£ c ) + (f 2 t?c + rizdr, -I- y 2 t? c ) 

(2.34) 

- + hiV’r? + (iV;) 

» (s2<£f + Vzd:, + + (fxXf + VxXr) + CzXc) 

( s y d -+ - \yi) ) 

The metric quantities used in the above transformations are defined as 
6- = ^ ( */r? - y;Zr,)- Vx = J(y s -2 t ‘ - yc2. : ) 

fy = J(x c 2 r/ Xy2,-), 7/y = J(xc2 s - - X<2c) 

f 2 = Jixrjy, - x-y v ), r) z = J(i c ye - x*y c ) (2.35) 

Crr = J{ycZri - Vr]Zc), Cy - J{x n Zc - XcZy) 

& = Jix^yr, - Xr,yc) 

where J is the Jacobian of the transformation given by 

J — Cxi^ySz ~ Vz^y) T hz(£zCy ~ fyfz) T Catfy^z — ^zVy) (2.36) 

Unsealed contravariant velocity components of the velocity vector q can be 
defined as 

U = £ X U + fyV + tzW 

V = y x u + T} y v + y 2 u> (2.37) 

W = $ x u + $ y v + $ z w 
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The contravariant velocity components are conveniently split between the 
potential contribution and the vector potential contribution as 

U = U 4 + U B 

v = v 4 + V B 
w = w 4 + W B 

where 

U 4 = A^<f>c + A^tr, + A 4< <f>< 

v 4, = A‘ v <j)£ + A^tr, + 

W 4 = A^4>c + A™<f>r, + A«<f> c 

U B = (V£ X Vf /) • Br, + (V£ X V?) • jB ; : 

V 5 = (Vtj x VO • Be + (Vtj x Vf) • B s - 

= (Vf x VO • Be + (Vf x Vg) • B v 

and 

4« = V£-V£ -£*+$ + £* 

4^ = Vg • Vr] = til -t r/j* + r\\ 

= Vf-Vf»fJ+ Sy -t S'? 

4^ = VOVr/ = i x r\ T + 4- 0r? 2 

4’’ ! ' = V£-Vf = £zfx + £j/fy 
4^ = Vr/ • Vf = T/jfx + ^j/ft/ + Wz$z 

It can also be shown that the operator g • V applied to any variable / 
expressed in the computational domain as 

(9 • V)/ = uf x + vf y + wfz 

= Uft + Vf r] + Wf c 


scalar 


( 2 . 38 ) 


( 2 . 39 ) 


( 2 . 40 ) 


can be 


( 2 . 41 ) 


Similarly, the operator Q • V which appears in the vorticity transport equation can 
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be transformed as 

(w • V)/ = f x + Ulfy + W3 f z 
— ^l/| + ^2 fr\ + ^3/c 

= (n- v)/ 

where 0 is the unsealed contravariant vorticity vector defined as 


Hi 

0 = I fi 2 
^3 


1 + iyW 2 + sc^’3 

r)xU\ -+ r}yU) 2 + VzU'i 
Cx^i + Cv w 2 + fz w 3 


(2.42) 


(2.43) 


L Governing equations in transformed coordinates 

The continuity equation can now be written in transformed coordinates as 

'p(A^4>c + A- r ><f>r, + a^ 4><) 


' p(A^<t>c + + AW<t> ( ) 


p{Ab<j>£ + A^<f> v + 


+ 


+ 


(2.44) 


U B \ j V B \ fw B \ 

7r s+ Tr ,+ 


The Poisson equations for the vector potential functions can be transformed simi- 
larly employing the same metric groupings and can be written in a compact vector 
form as 


( A^Bc + A^Br, + A&B S 

V " 7 " 

/ A^Bc + A^Br, + A™B f 

V ' 1 


+ 


+ 


(2.45) 


A^Bc + A 7 * Bn + A^BA _ w 

J / ~~ ~ J 
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In the same transformed coordinates the Crocco relations are given as 
vw 3 - wu>2 = - T(£ X S£ + rixSf) + $ X S{) 

U’Ul - uw 3 = - T(tySc 4 TJySrj 4 Cy5 s *) 

uu.’2 ~ = - T(£ z sc 4 r) 2 Sfj + 

while the vorticity transport equation (2.32) is written as 


p{UuJc + Vu n = — 


' A^Qc 4 A^Ur, 4 A^ 


Ur 


A^uc + A^Ur, 4 
- 

A^Uc 4 A^Un + A^lUr ^ 


J 


4 


(2.46) 


4 (H V)<f 

(2.47) 


and the entropy convection equation is given as 


U sc 4 V Sy 4 W — 0 (2.48) 

The consistency equation for the vector potential functions takes the form 

V • B = • Be + Vri • B n 4 Vc • B s - = 0 (2.49) 

or simply 

(^z^c + Vx^T] + + 

(CyXf + VyXy + SyX<) + (2.50) 

(fz^f + + $*^f) = 0 

Similarly, the vorticity consistency equation (2.18) can be expressed as 

[fx( w l)| + Vxiuijy 4 £c(wi)<] + 

Uy( w 2)£ + J 7y( u; 2)t) + Cy( w 2k] + ( 2 -5l) 

+ £s( w 3)?] = 0 


E. Numerical Algorithms 

The dual potential formulation outlined in the previous sections has several 
advantages. In an iterative solution scheme the governing equations (2.44) - (2.51) 
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are weakly coupled and hence can be solved separately. Because of the uncoupling, 
efficient solution schemes can be used for solving each individual equation. The 
actual iterative solution strategy for inviscid and viscous flow will be presented in 
Chapters III and IV. The numerical algorithms for the individual equations are 
described below. 

The Poisson equations and the continuity (full potential) equations are con- 
ventionally differenced (i.e., central) and are solved with ADI-like procedures that 
use a sequence of relaxation parameters. The convection equations are differenced 
using upwind differencing in x or f and central differencing in the other directions. 
An implicit 2-factor approximate factorization (AF) method is being used as the 
relaxation scheme. Tangency and boundary consistency relations are enforced using 
one-sided differencing in the direction normal to the boundary surface and central 
differencing in the other two directions on the boundary surface. An ADI procedure 
is also used to relax the consistency and tangency boundary values. The following 
notation is used in describing the numerical algorithms. In the uniform computa- 
tional domain, Af = Ar? = Af = 1. The indices j,k,l refer to grid points in the 
£-, ??-, and f-directions, respectively. Generally a subscript is not shown unless it 
varies; for example, The difference operators are defined in terms 

of the shift operator Ef ] Uj = Uj±j as 

V ? =(l-£7>)/(Ae) 

A { = (Et 1 - 1)/(AJ) 

= (E+ 1 - E7')/(2A() 

« { = (E^-E'b/( AC) 


(2.52) 
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Second-order accurate one-sided difference operators are defined as 



(2.53) 


T. ADI algorithm for the continuity and Poisson equati ons 

The continuity equation for 4>, equation (2.44), is updated using an approximate 
factorization (ADI-type) algorithm in delta-form. The steady-state form of the one 
described in Bridgeman et al. (1982) is used and details can be found there. The 
differencing scheme, in AF form, is given by 

i / - hbc{pA^) n lc]\I - h6 ri (pA T,r, ) n 6 r) )\I - h'6,(pA^) n 6.\ 

> (©"" 1 4> n ) = hu [*c(pf/*) n + f>r,{pV 4 V 4- k{j>w*) n (2.54) 

U B V B c W B c _ 1 

+ —J- bcp 4 ~j- Orfp 4 -J—O S -P - /too 

where p — p/ J , and the superscript n refers to the iteration level. The differencing 
of the contravariant velocity component terms that appear on the right-hand side 
of equation (2.54) is illustrated by considering and U B . Central differencing in 
space is used throughout, except for 4> terms which are associated with a second 
derivative in f. Thus, 

A e *, + 4 (a%M^ + 

+ \ 1 + A ?W,) 

and 

u B = i x u B + i y v B + i z w B 


(2.55) 


(2.56) 
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where 

u B = + V ySr,^ + $ y 6 C0) - (6^X + rizSrjX + f z 6 ( \ ) 

v B = {Sz&cd + r] z 6 v d + c *M) - (txtcip + rj x 6r,rp + $ x 6$il>) (2-57) 

W B = (Zx 6 £X + rj x 6riX +■ Cx^-x) - UyAtti + »7jA,i? + ^6-tf) 

The V and IV terms are also treated with central differencing except for (f> terms 

associated with second derivatives in /? alone and f alone, respectively. 

The free-stream subtraction term. /?oc- appearing in the right-hand side of equa- 
tion (2.54), accounts for incomplete metric cancellation (Bridgeman et al., 1982) and 
is given by 


r, l ( Po oU<x\ , / Poe f'oo 

k°° = h {-j- ) + 1”, {- — 


6 C 


^ PooW, 


(2.58) 


An ADI algorithm is used to solve equation (2.54). Rewriting equation (2.54) 
in the form 

LcLr,L<(<t>”- n - 0 n ) = R (2.59) 


the ADI algorithm is implemented in three steps as 

LcA(f>* = R 

L n A<t>** = A (f>* 

L'A<f> n+1 = A4>** 


(2.60) 


4> n+1 .= </> n + Acf> n 

The algorithm given by equation (2.59) requires only a series of scalar, tridiago- 
nal inversions and is, therefore, solved efficiently. It may also be noted that the 
rotational terms only appear explicitly in the residual. The relaxation parameter 
h, which appears in equation (2.54), is determined from the geometric sequence 
(Ballhaus et ah, 1978): 


A = 



1 


, i = 1,2,3... N 


(2.61) 
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where A = 1 Jh. A second overrelaxation parameter u> is used to scale the residual. 

Each of the Poisson equations for the vector potential functions is updated in a 
similar manner, which results in the AF form, 


/ - Me 


x(£ n+1 - B n ) = huj 


. 1 1 
A'''* 
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- ( xm \ n _ 1 r .. / A''> \ n - 1 
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where 
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+ 5 KiM'+i + Afi ; B,) 


- 4^ ! A,,B t 4 - (4£+,^Bfr + , 4 


B 


l”* A) 

+ 5 (4? +1 *,B, + , + 4? J «,B,) 


(2.62) 


(2.63a) 


(2.63b) 


The scheme given by equation (2.62) is implemented by the same three-step ADI 
algorithm given by equation (2.59). 


2. Convection differencing and relaxation 

The differencing used with the convection equation is illustrated by considering 
the entropy convection equation in general coordinates: 

Us e + Vs n + Ws^ = 0 (2.64) 

Assuming U is generally larger than V and W , this first order wave equation is 
discretized by using three-point, second-order-accurate central differencing in the 
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r?- and ^-directions, and three-point, second-order-accurate upwind differencing in 
the ^-direction. Upwind differencing of Use is automatically achieved using 

Use = U+6 h c s + U~6 f c s (2.65) 

s ^ S 


where 


17+ = M and U- = U -' U ' 


2 2 
The convection equation is then solved using the AF relaxation algorithm 

(/ + hU + 6 h c + hW6 s -) (/ + hU~S i i + hVS^j (s n+1 - s n ) 


( 2 . 66 ) 


-/i (u + 6 b c + U~bl + V6 v + W«.) s n 


(2.67) 


where h is another relaxation parameter ( h > 0). Because central differencing 
is used in the r?- and ^-directions, it is necessary to add fourth-order numerical 
dissipation terms to the differencing scheme to control \he odd - even uncoupling 
of grid points and to control any nonlinear effects such as shocks (Pulliam, 1985). 
This is done explicitly in the numerical algorithm in order to be able to invert only 
tridiagonal matrices in the rj- and ^-directions. Second-order implicit dissipation 
is added to stabilize the explicit fourth-order dissipation for large values of the 
relaxation parameter h. Adding the numerical dissipation terms to the differencing 
scheme given by equation (2.67) results in the following scheme for the convection 
equation. 

(/ + hU+6% + hWS- - 3h|VF|AV| f )(7 + hU~ 6^ + hV - 3/r|U|AV|,,) 

x(s n+1 - s n ) = -h [U + 6 b c + U~6* + V6r, + VF6< + |V|(AV) 2 |^ + |W|(AV)\]s n 

( 2 . 68 ) 

Since the second-order dissipation terms are added implicitly, they do not change the 
steady-state solution of the convection equation. The algorithm given by equation 
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(2.68) is implemented in two steps as 

(/ + hU + 6^ + h\V6 < - 3/ijiy|AV| ? )As* = 

' (2.69) 

-h [U + S h c + U~6l + V6 n + W6, : + jVj(AV) 2 !^ + ;W|(AV) 2 !,]s n 

(/ 4 hU~s! -f hVdrj - 3fcjV:AV| r/ )(s n+1 - s n ) = As* (2.70) 

In the first step, a forward sweep in £ is achieved by requiring a series of tridiagonal 
inversions in f. In the second step, a backward sweep in £ is achieved by requiring 
a series of tridiagonal inversions in 77. 

A similar algorithm can be used for solving the vorticity transport equations 
(2.47) by replacing the numerical dissipation terms with the viscous diffusion terms. 
The details of the algorithm used for solving the vorticity transport equation will 
be presented in Chapter IV. 


3. Integration algorithm for vorticity 

The consistency relation for vorticity, equation (2.51), is used to integrate u.^ 
in the ^-direction from an upstream boundary to the outflow boundary. Equation 
(2.51) can be represented in the form 


aifc + a 2 fr, + 03/5: — a 4 


(2.71) 


where / = u/j, aj = £*, «2 = a 3 = £z, and 04 contains the remaining terms of 
equation (2.51). 

Consider differencing this equation on a j-plane using first-order one-sided dif- 
ferencing in £ and central differencing in r? and f. That is, at the j th plane 

(/fc+l ~ fk- 1) | 

3 (2.72) 


°1 {fj ~ fj- 1) + a 2 


(fl+\ ~ fl-l) > 

+ «3 ~ \j = a 4 
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Putting into the delta form, fj - fj-\, dividing by aj, and approximately factoring, 
this difference equation is written as 


u + 


~ fj- 1 ) - 

O) a\ 

a 4 ,«2r a 3,w 

—A, - - f > :)/,_! 

a, a j fl] J 


(2.73) 


If three-point, second-order forward differencing is used for /^, this becomes 


2 a j 


: a ; 


(' + 3a,W -sfWi -/,■-.> = 

2 G4 1 . 2 do 

3^ 4 A- 1 " / '- 2 > _ s ( «A + 1 


(2.74) 


Note that this approximate factorization is using a space delta, not an iteration 
delta. First-order backward differencing is used at the first £ = constant plane, and 
second-order differencing is used at all successive planes to integrate the consistency 
relation for uq . At each £-plane, the algorithm requires a series of simple tridiagonal 
inversions in ?] and f. 


4_. Consistency boundary condition 


The consistency relation for the vector potential functions is implemented at 
boundary surfaces by a procedure similar to the one described above for integrating 
the vorticity consistency relation. Consistency boundary conditions are illustrated 
below for the case in which the body boundary-surface coincides with af = constant 
plane. Other boundary surfaces receive similar treatment. 

Differencing equation (2.50) on the £ = constant surface, 1=1, using first-order 
one-sided differencing in and central differencing in £ and rj, gives 


Zz ^ !/=i 

+ t \i = i + fz(t/>2 - 0l) = «4 i/=l 


(2.75) 
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(2.76) 


where 

a 4 = + fxA^t?) 

+ (M fX + ri y 6 n x + CyAsX)] 

Putting this into the delta form, dividing by — and approximately 

factoring, the difference equation 


c 

S .2 


[l - ^6;)(I - !!%)(*, - *) = + ^S v U-2 - ^ (2.77) 

Sz Sz SZ S z S z 


is solved on the surface for (0] — fa) by inverting tridiagonal matrices in f and 
then in r). If three-point second-order forward differencing is used for c-derivatives, 
equation (2.77) is modified in the same way that equation (2.73) is changed to 
equation (2.74). 


5. Tangency bo undary c ondit ion 

Tangencv or no-flow-through is imposed on a boundary surface by setting the 
appropriate contravariant velocity component to zero. On a ( = constant plane, 
W is set to zero. Tangency is enforced through implicit boundary conditions on <f > , 
which are obtained by solving the continuity equation at the half-cells neighboring 
a solid boundary. The procedure is illustrated by considering a finite-difference cell 
(Figure 1). The center of the cell is located at ( j,k,l + |). The finite-difference 
form of equation (2.44) can be written as 

+ 

+ 

= 0 


W U >j+k,kMi~ 

A/7 

(Wj, <:,/+! 

AC/2 


(2.78) 




Figure 1. Finite-difference cell near a wall boundary 
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Evaluating variables at l + \ by weighted averaging between l and / + 1 as, for 


example, 


3 1 

= (i + 


(2.79) 


where 


v i+> ,u =■ 

1 Cs c* 

+ -(j4yj, j A.*d> J+ j + AjAc<p]) f vector potential terms 


(2.80) 


equation (2.78) can be written as 


(i + ^ f )(pc) J+r (/+ 

(I 4 iA f )(^) t+ , - (/ 4 t A,)(pV) fc _ , + 2 (pW) l+L = 0 


(2.81) 


To facilitate the application of approximate factorization, the cross-derivative and 
vector potential terms are lagged in time in the usual way to obtain the relaxation 


algorithm 


|/ - - hi n WTK\\ i - , A f 

x (<t, M 1 -</,-) = + i„(pV) n 4 2 {p\vy^, 


Equation (2.82) is of the same form as equation (2.54) and hence tangency can be 
enforced implicitly in the ADI algorithm given by equation (2.59) with the tridiag- 
onal and right-hand-side terms modified appropriately. 
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III. INVISCID ROTATIONAL FLOW 

For many practical flow situations, the inviscid approximation gives a very useful 
and accurate estimate of the flowfield. This is especially true when viscous effects 
are small and confined to a thin shear layer adjacent to the wall and when the flow 
remains attached and does not separate, thereby diffusing and convecting vorticity 
into the outer flow. In this chapter, the solution scheme for the Euler equations is 
presented together with applications to inviscid three-dimensional rotational flow. 
The approach taken is to increase the complexity of the approximations regarding 
the flow. First, the solution to an irrotational inviscid flow is obtained. Then, 
an inviscid solution of the secondary flow produced by a nonuniform stagnation 
pressure is obtained. In particular, the three-dimensional dual potential procedure 
outlined in Chapter 11 has been used to simulate the flow through the inlet of an 
indraft wind tunnel and numerical results for this case are presented in this chapter 
along with some experimental data for comparison. 

A. Solution Strategy 

The dual potential formulation outlined in the previous chapter has several 
advantages. In an iterative solution scheme, the governing equations (2.44) - (2.46), 
(2.48) - (2.51) are weakly coupled and hence can be solved separately as described 
below. The rotational and irrotational components of velocity are also decoupled. 
As a result, the vector potential functions need be determined in only the rotational 
part of the flow domain. 

Because the equations are weakly coupled, an efficient iterative solution pro- 
cedure can be used with each equation updated sequentially using the numerical 
algorithms described in Chapter II. At each grid point an initial guess is made for 

HiECEDINQ PAGE UMfd W#T KIK'-O 
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values of 0, ti, 0, s, p, and for iuj, tu 2 , and u> 3 . The following iteration scheme is 
then implemented: 

1 ) Vector potential functions, t?,y,0 . Individually update tf,y, and 0 from the 
Poisson equations given by equation (2.45) for assumed values of uq,w 2 , and W 3 . 
Boundary conditions for these equations are kept compatible with the consistency 
relation, equation (2.50). 

2 ) Scalar potential funct ion, 0. Update 0 from the continuity equation (2.44) 
using the ADI-like algorithm for the transonic full potential function and previously 
updated values of tf,x, and 0 . 

3) Entropy, s . Using updated values of $, y, 0, and 0, evaluate U, V , and VP and 
update s from the convection equation (2.48). 

4) Den sity, p . Update density from the Bernoulli equation (2.28), using previ- 
ously updated values of u.v,w , and s. 

5) Vorticities, u.'i,u; 2 , and u$ . Vorticity components are evaluated from Crocco’s 
equations and the consistency relation for vorticity. Assuming u > |v| and |ip| and 
u ^ 0, w 2 and W 3 are evaluated from equation (2.46) as 

U > 2 = [VU\ ~ T{£ z Sc + TJ z Sr, + feS,-)] ! U 

(3.1) 

W 3 = [tuwi + T{i y sc 4- rjySy ■+ ; y s . : ) J ju 

and u>\ is determined from the vorticity consistency relation, equation (2.51). 

6 ) Test for convergence and, if necessary, return to (l). 

For incompressible flow, density is taken to be constant, and Crocco’s equations 
and the entropy convection equation are replaced by Lamb’s equations and the total 
pressure convection equation, respectively. In this case, u ; 2 and W 3 are evaluated as 
w 2 = [ vuJi + UzPO' + rj z p 0rj + $zPQf)/p]/u 

w 3 = [u)W] - Uj/Pq^ + r) y p 0r] + $yPQ.)/p}/u 


(3.2) 
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and u>i is again determined from the vorticity consistency relation. 

B. Geometry and Grid 

As a first application, the three-dimensional dual potential procedure outlined 
above was used to compute the flow through the 80- by 120-ft leg of the National 
Full Scale Aerodynamics Complex (NFAC) at the NASA Ames Research Center. 
A view of the complex is shown in Figure 2. Screens and vanes are placed at the 
entrance of the inlet in order to isolate the test section from outside winds. Also, 
the screens prevent birds and other objects from entering the tunnel. The presence 
of these screens and vanes is accounted for in the numerical simulation as a jump 
condition based on actuator disk theory. The contraction ratio for the tunnel is 5:1 
and the Mach number in the test section does not exceed 0.2. The fluid is, therefore, 
assumed to be incompressible in the numerical procedure. 

A three-dimensional grid to model the 80- by 120-ft leg of the NFAC was gen- 
erated using simple algebraic shearing. Figure 3 shows a view of the tunnel walls 
and floor. A comparison of Figures 2 and 3 shows that the grid represents the 
NFAC accurately except for the rounded lip at the entrance of the inlet. To avoid 
a complex grid generation task, the rounded lip was not modeled. Since the 80- 
by 120-ft wind tunnel rests on the ground, the grid and the boundary conditions 
must model the ground plane in order to obtain the correct inflow conditions. The 
computational domain includes the wind tunnel and its surrounding region to facil- 
itate the application of boundary conditions in the far-field. All calculations were 
performed on a 39 x 41 x 34 grid of which 26 x 17 x 21 grid points are interior to 
the wind tunnel. The tunnel walls are modeled to have a thickness equal to a single 
grid increment. The actuator disk modeling the screens and/or vanes is located at 
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Figure 3. 


Perspective view of the grid for the NFAC inlet model showing 
the ground plane and tunnel walls 
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a f = constant plane corresponding to the location of the trailing edge of the vanes 
in the actual configuration. Figures 4 and 5 show views of three different planes 
of the grid. A detail of the grid is given in Figure 5 to illustrate the modeling of 
the tunnel walls. The average width of the test section, 100 feet, was chosen as the 
reference length in the calculations, and all results are presented accordingly. 

The flow is assumed to be irrotational outside of the tunnel, and only the 
potential function, <£, needs to be computed there. Upstream and in t he farfield, the 
value of the potential is kept constant at its free-stream value. At outflow, outside 
the tunnel, the streamwise component of velocity is allowed to vary by using a 
zero-gradient boundary condition on U . Tangency is enforced at the ground plane 
and at all the tunnel walls (interior and exterior) through the implicit boundary 
conditions on 4> described earlier. The flow rate through the tunnel is determined by 
specifying the axial velocity <p x at the outflow boundary within the tunnel. On each 
tunnel wall, the vector potential corresponding to the near-normal direction to the 
surface is obtained from its consistency relation, and the other two vector potential 
functions are kept constant. For example, on the ground plane f = constant, d and 
X are set to zero, while ip is obtained from equation (2.50). At the outflow boundary 
£ = constant, \ and ip are obtained using a second-order extrapolation technique, 
whereas d is computed from the vector potential consistency relation. 

In the absence of any screens and/or vanes, the entire flow is assumed irrota- 
tional and hence only <f> needs to computed, with the other variables being constant. 
With screens and vanes, a total pressure loss across these devices sets up a rota- 
tional flow field inside the tunnel which is accounted for by solving for the total 
pressure (incompressible formulation), vorticity, and the vector potential functions. 
It is necessary to solve for these variables only within the tunnel since the flow 
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(a) 



(b) 


Figure 4. Cross-sectional views of the sheared grid for the 80- by 120-ft leg 
of the NFAC. a) Grid in the y-z plane at inlet entrance; b) Grid 
in the x-z plane at mid-span 
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Figure 5. Cross-sectional view of the sheared grid for the 80- by 120-ft leg 
of the NFAC in the x-y plane at mid-height 
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upstream of the vanes and screens is taken to be irrotational. Inflow boundary 
conditions for total pressure and the vector potential functions are obtained from 

the actuator disk model for the vanes and screens. 

I 

| C. Actuator Disk Model 

| The presence of vanes and screens at the entrance of the wind-tunnel inlet is 

modeled using actuator disk theory (Horlock, 1978; Ross et al., 1986) in which these 
devices are idealized to a jump condition. Vanes and screens constrain the flow in 
some specified direction and produce a drag force opposing the motion of the fluid. 
The drag or resistance produced by vanes and screens can be expressed in terms of 
the total pressure loss across the actuator disk model. For a given screen or vane, 
the total pressure drop can be determined empirically as a function of the local 
dynamic head as 

P0 2 = POi - K \pqi ( 3 - 3 ) 

Here, the subscripts 1 and 2 refer to locations immediately upstream and down- 
stream of the actuator disk as indicated in Figure 6. The loss coefficient K is 
empirically determined as a function of the local Reynolds number, onset flow an- 
gle, and the geometrical characteristics of the screens and vanes. The contribution 
to the loss coefficient from the vanes and the screens may be separated into two 
components as 

K = K v + K s (3.4) 

where K v denotes the vane-loss coefficient and K s the screen-loss coefficient. 



ACTUATOR 

DISK 


Figure 



P0 2 = P0 1 -K s -lp(u2 + v 2) 


v 2 = (1 -C) Vl 

K s & C = func(Re, a^) 


6. Deflection of flow through an actuator disk model for a screen 
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T Flow through a screen 

The screen is modeled as an actuator disk located in the flow field as indicated 
schematically in Figure 6. Variables at the locations upstream and downstream 
of the disk are denoted with the subscripts 1 and 2. respectively. The following 
equations describe incompressible flow through the disk. 

Continuity: 

U] = u-2 (3.5) 

Momentum in x-direction: 

Pl - P 2 ~ D x = pul ~ pu\ (3.6) 


where D x is the drag force in the x-direction and is obtained from the characteristics 
of the screen and the upstream flow conditions. The loss coefficient, based on the 
static pressure drop and expressed as a function of the dynamic head is found to 
be a function of the local Reynolds number, (Re), for onset angles, a, up to 40°. 
That is, 


K 


D x 

$ ° ~ 1 2 

2 PT 


Pl ~ Pl 

\pq 2 


= f(Re) 


(3.7) 


In addition to creating a pressure drop, the screens have the effect of straight- 
ening the flow through them. It is found (Horlock, 1978) that at a given Reynolds 
number, the change in the tangent of the flow angle is a function of the upstream 
flow angle, that is, 


tan ai — tan a 2 = /(tan ax) 


(3.8) 


or 


u\ tan ax — «2 tan a 2 — F{ui tan ax) 


(3.9) 
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For onset flow angles up to 40°, the change in the tangent of the flow angle is found 
to be proportional to the upstream flow angle, that is 

c\ — C 2 — Cc\ (3.10) 

where c - v tan a = (v 2 + u ,2 ) 1/2 is the lateral component of velocity. 

The empirical parameters, K so and C are determined from the upstream con- 
ditions and the physical screen characteristics (porosity, w'ire diameter). An addi- 
tional simplification is achieved by assuming that equation (3.10) is valid for both 
tangential components of velocity, that is, 


v l ~ v 2 = C\)\ (3.11) 

uq - w 2 = Cu>i (3.12) 

The total pressure drop across the screen is given by 


which reduces to 


POj - Po. 


1 

2 


pq 2 


P\ Pi 

\ p <! 2 



K< = K S o + (2 C - C 2 ) sin 2 aj 


(3.13) 


(3.14) 


The two constants, K$ 0 and C needed to determine the total pressure loss for 
the NFAC have been determined experimentally and reported by Ross et al. (1986) 
and van Aken (1986). 


2^ Flow through vanes 

In addition to screens, many wind tunnels employ a a cascade of vanes at the 
entrance in order to assist in providing a uniform test-section flow. The distribution 
of the splay angle of the cascade is an important factor in minimizing the variation 
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in test-section total pressure (Ross et al., 1986). The loss coefficient for the cascade 
is a function of the individual drag coefficient of the vanes including its geometric 
characteristics such as the pitch, chord and solidity. For the 80- by 120-ft leg of the 
NFAC, the loss coefficient for the cascade is obtained from van Aken (1986) and is 
found to be 0.4. Another important effect of the cascade is to direct the flow at 
the inlet entrance at a specified angle thereby determining the lateral component 
of velocity as 

V'2 = U'2 f an 0 (3.15) 

where 6 is the splay or turning angle. Various splay angle distributions were tried 
for the NFAC and have been reported by Ross et al. (1986). Finally, it may 
also be useful and even necessary, as in the case of the NFAC to control the vertical 
component of velocity, w. This is done by the construction of splitter plates, running 
laterally between the two side walls of the tunnel, which are located at different 
heights of the vanes to form something like a honeycomb screen at the inlet entrance. 
For the NFAC the splitter plates have been constructed horizontally along the width 
of the tunnel, so that the fluid leaves the plates without any vertical component of 
velocity, that is, 

w 2 = 0 (3.16) 

The coefficients K v , K?, and C used in this study are based on the experimental 
data of Ross et al. (1986) and van Aken (1986). Except for one case (Figure 9), 
values of K s = 1.8 ,K V = 0.4, and C = 0.2 were used in the numerical simulations. 

3. Implementation of actuator disk model 

This actuator disk model is implemented in the numerical simulation by adopt- 
ing the following procedure at each iteration level. Letting j be the grid point at 
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which the actuator disk is located, and using the subscripts 1 and 2 for values just 
upstream and downstream of the disk, 

• at j are calculated from equation (2.34) using backward differences 
for ihe x- or ^-derivatives, that is. making use of values at j and j - 1 . 

• po al J ' s calculated from equation (3.3) using appropriate values of K . It may 
be noted that po 1 is a constant for the irrotational flow upstream and hence 
known from free-stream conditions. 

• V 2 and u >2 at j are calculated using equations (3.11) and (3.12) for screens alone 
or equations (3.15) and (3.16) for screens and vanes. 

• toi at j is calculated from its definition using values of V 2 and W 2 as 

— Wy - v z (317) 


d is set to zero at j. 

X and at j are calculated so that i >2 and are prescribed using forward 
differences in x as 


= V’j-H - [4>y - V 2 )&x 
X] = Xj + 1 - [w 2 - 4>z)Ax 


(3.18) 


• Finally, the continuity equation through the actuator disk, equation (3.5), is 
enforced by evaluating 4> at j using appropriate differencing in x. That is, 
equation (3.5) is differenced as 


“1 = {<t>j ~ <f>j-i)/Ax = u 2 = - <t>j)l Ax -j- {xl> y - Xz) (3.19) 

Collecting terms involving 4> on one side, the above equation can be rewritten 
as 


( - ^>+ 1 + 2( t>] ~ 4>j- 1 ) = {4>y - Xz)Ax 


(3.20) 
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This constraint on <f> 3 is enforced implicitly in the ADI algorithm for 4> to ac- 
celerate convergence and improve stability. 

D. Results and Discussion 


T lrrota t ional flow 

As a simple verification of the numerical procedure, Figures 7 and 8 show so- 
lutions obtained for potential flow alone (K = 0) compared with other numerical 
computations and experimental data as taken from Kaul et al. (1985). The exper- 
imental results presented throughout were obtained from a 1/ 15-scale model of the 
NFAC. Figure 7 shows the variation of the pressure coefficient, Cp = (p-p^/^puf , 
on the side wall versus length at mid-height of the tunnel. Here, the subscript t 
refers to the test section, which is at the outflow boundary in the present simula- 
tion. Experimental data are also plotted in Figure 7 and the agreement between 
the present results and the experimentally measured values is quite good. The spike 
in the present results is likely a result of the simple lip treatment. Figure 8 is a 
comparison of the pressure coefficient along the center of the tunnel floor, with Eu- 
ler and panel methods. The panel method results agree very well with the present 
results, whereas the Euler method differs in the entrance region of the inlet. These 
calculations required about 5 //sec per iteration per grid point on a CRAY XMP 
vector processor and fully converged results were obtained in 900 iterations which 
correspond to about 240 sec of CPU time. 

2. Rotational flow 

In the design of wind tunnels, it is, of course, essential to maintain uniform 
flow in the test section, and screens and vanes (located upstream of the test sec- 
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Figure 7. Variation of pressure coefficient with length along a side wall at 
mid-height of the tunnel 
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Figure 8. Variation of pressure coefficient with length on the tunnel floor 
at mid-span 
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tion) are frequently used to improve the flow quality in the test section and isolate 
the wind tunnel from external turbulence. However, these devices induce a total 
pressure loss that varies with the upstream dynamic pressure. If they are placed in 
a region with substantial variation of upstream velocity, they can induce a signif- 
icant total pressure variation (and, therefore, rotational flow) that can persist far 
downstream. Placing a screen in the wind-tunnel inlet can therefore induce an un- 
desirable nonuniform flow in the test section. This effect is demonstrated in Figure 
9 for the original NFAC inlet; the figure shows the variation of the dynamic pressure 
in the test section with the spanwise distance, y, at mid-height in the tunnel. The 
quantity plotted in this figure is the normalized dynamic pressure, (g/<7 c /) 2 , where 
the subscript cl refers to the centerline. Three numerical results are shown in the 
figure: screens alone, screens and straight vanes, and screens and splayed vanes. In 
one case, the vanes are straight, that is, they are aligned with the axis of the tunnel 
\6 = 0 in equation (3.15)] just as in the original NFAC design. In the other case 
the vanes are splayed at angles developed for the redesigned NFAC (Ross et al., 
1986; van Aken, 1986). Along with the numerical results, experimental data from 
a 1/15-scale model with straight vanes (van Aken, 1986) are also plotted in this 
figure. The agreement with experiment is relatively good and indicates that the 
original inlet configuration results in an unacceptably high variation of flow quality 
in the test section. 

A numerical and experimental study of the NFAC (Ross et ah, 1986) resulted 
in a slight modification to the inlet geometry. The curved walls in the span (or 
y-direction) were replaced with straight walls near the entrance of the inlet and 
the vanes were splayed at angles that are close to the “ideal” streamline angles 
at the entrance of the inlet. A schematic view of this new geometry is shown in 
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Figure 10; again, the true lip geometry is not represented. Results obtained for 
this modified tunnel are shown in Figures 10-14. Also shown in these figures are 
experimental data from a 1/15-scale model (van Aken, 1986). Figure 10 shows the 
variation of the total pressure coefficient, (po"Po /)/ \pQdi with spanwise distance at 
mid-height of the test section. Again, two splay distributions were examined: zero 
splay angles and the redesigned splay distribution. The results show a dramatic 
improvement in test-section flow quality for the case of splay distribution and the 
agreement with experiment is good except near the walls. This is perhaps expected 
because viscous effects are not taken into account in the present procedure and the 
correct lip geometry is not modeled. Figures 11 through 14 compare predictions 
with measurements of the dynamic pressure at various height and span locations 
in the test section. Overall, the agreement with experiment is satisfactory except 
near the walls. 

The above calculations were performed with zero free-stream velocity, that is, 
with no wind outside the tunnel. With the present code, a variety of outside wind 
conditions can be simulated by changing the far-field boundary conditions on 4> (in- 
cluding rotational incoming flow if the vector potential functions are solved there). 
The effect of wind direction on the test-section flow quality was studied by simulat- 
ing a case in which the wind velocity was about 15% of the test-section velocity and 
blowing at an angle of 45° with respect to the axis of the tunnel. The results are 
presented in Figure 15, along with results for the case in which the wind was blowing 
at 0°. Figure 15a shows the spanwise variation of dynamic pressure at mid-height, 
and Figure 15b shows the variation of the dynamic pressure with height at mid-span 
of the test section. The results indicate that the vanes, inlet length, and screens 
do a relatively good job of isolating the test section from such extreme variation in 
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Figure 15. Effect of wind direction on dynamic pressure in the test section. 

a) Spanwise variation of dynamic pressure at mid-height; b) Vari- 
ation of dynamic pressure along height at mid-span 
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external flow conditions. This has been confirmed experimentally though quantita- 
tive data are not available for comparison. Figure 16 shows the velocity vectors in 
the horizontal and vertical planes at midspan and midheight, respectively, for the 
case in which the wind is at an angle of 45°. 



Figure 16. Plot of velocity vectors for wind blowing at 45°. a) Velocity vec- 
tors at mid-height; b) Velocity vectors at mid-span 
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IV. EVALUATION OF A DUAL POTENTIAL - 
BOUNDARY LAYER INTERACTION SCHEME 

The dual potential decomposition scheme was applied to compute inviscid irro- 
tational and rotational flows in Chapter III. It is possible to use the same formulation 
to compute viscous flows by additionally solving equation (2.47) for the transport of 
vorticity. Most numerical techniques for equations of this type are notoriously inef- 
ficient because of the multiple length scale phenomena that arise in high Rey nolds 
number viscous flows. With the dual potential formulation presented earlier, it 
might be possible to take advantage of the speed and efficiency of boundary-layer 
algorithms to compute the vorticity - which can be used to determine the veloc- 
ity field by solving the Poisson equations for the vector potential functions. The 
validity of such an interaction scheme for separated flow is examined in a prelim- 
inary way by computing the flow over a two-dimensional configuration which has 
been studied by others using conventional viscous-inviscid interaction schemes. The 
present dual potential interaction scheme is also applied to compute the attached 
flow over a three-dimensional trough configuration. These three-dimensional results 
are compared with other numerical data. 

The three-dimensional boundary-layer equations are presented in §IV.A. A re- 
duced version of the dual potential formulation for two-dimensional viscous flows is 
also presented in this section for convenience. The interaction scheme between the 
dual potential and boundary-layer equations is discussed in §IV.B. Applications of 
the present interaction scheme to two- and three-dimensional trough configurations 
are presented in §IV.C and §IV.D respectively. 

All the calculations presented in this chapter were made with the flow assumed 
to be incompressible (density was assumed to be constant). 
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A. Governing Equations for Viscous Flow 
T T hree-dimensional boun d ary-layer e qu atio ns 

The boundary-layer equations are derived from the Navier-Stokes equations 
using an order of magnitude analysis (Schlicting, 1979). The basic assumption in 
deriving the boundary-layer equations is Prandtl’s hypothesis that for sufficiently 
large Reynolds numbers, there is a thin layer adjacent to the wall where viscous 
effects are at least as important as inertial effects. It is found from an order of 
magnitude analysis that viscous diffusion terms in the streamwise direction are 
negligible. Further, it is found that pressure variations in the normal direction are 
negligible and the normal momentum equation need not be solved. To distinguish 
the boundary-layer variables from the dual potential ones, boundary-layer variables 
are indicated with bars (“). In a regular Cartesian coordinate system, the boundary 
layer equations for steady, incompressible laminar flow are (Anderson et al., 1984): 
continuity: 

Ux + Vy + Wz = 0 (4.1) 

streamwise, i-momentum: 

p(uux + vug + wiiz) = -px + iizz (4.2) 

cross-flow, y-momentum: 

P{UVX + VVy + WVz) = ~py + Vzz (4.3) 

In these equations, the variables are nondimensionalized by free stream values as 
indicated in Chapter II. Also, the normal coordinate z and the corresponding veloc- 
ity component w are scaled by the square root of the free-stream Reynolds number 





Re so that they may have the same magnitude as the streamwise components (for 
laminar flow). 

a. Boundary conditions At a solid wall boundary, the no-slip boundary con- 
dition was applied: 

v(x,y,0) = v(x,y,0) = u>(x,y,0) = 0 (4.4) 

If no streamwise separat ion is present, these equations can be solved with pres- 
sure specified as the forcing function (the direct mode), where the pressure is cal- 
culated from the inviscid flow variables. Near and in reversed-flow regions, the 
boundary-layer equations are solved in the inverse mode to avoid saddle-point be- 
havior at the separation point. To solve the momentum equations in the inverse 
mode, the boundary conditions are modified so that inverse forcing functions can be 
specified. In conventional viscous-inviscid interaction approaches, the inviscid flow' 
is computed over a equivalent body rather than the actual body surface, that is, 
viscous effects are assumed to have the effect of altering the effective body surface 
by an amount equal to the displacement thickness and it is perhaps convenient to 
specify the displacement thickness as the inverse forcing function ( cf Carter and 
Wornum, 1975; Kwon and Pletcher, 1978; Lee and Pletcher, 1986). However, fol- 
lowing the work of Klineberg and Steger (1974), and Van Dalsem and Steger (1983), 
the wall shear stress is specified as the inverse forcing function in this study. The 
inverse forcing function is specified in the solution scheme by replacing the direct 
forcing functions (pressure terms) with expressions containing the inverse forcing 
functions. In the present procedure, these relations were obtained by evaluating 
the momentum equations at the wall. For example, the x-momentum equation 
evaluated at the wall gives 


Px — u zz 


(4.5) 


A Taylor series expansion about a point at the wall gives 

u 2 = «1 + U £ (z2 - ^l) + ~Uzz{z2 - ^l) 2 + higher order terms (4.6) 
Defining the nondimensional shear stress at the wall as 


T U)X 


u 2 


U’ 


(4.7) 


and making use of equations (4.5) - (4.6) gives the following relationship between 
the pressure gradient and the component of wall shear stress in the x-direction: 


Px — 


2 


U2-U1 r 
Z 2 -Z l Twx 


*2 - Z\ 


(4.8) 


In this work, the wall shear stress (r UJX ) was used as the inverse boundary 
condition, and the inverse method was used in only the x-momentum equation, 
with the y-momentum equation solved in the direct mode. The inverse method can 
be used even for attached flows, though it was used mainly in the separated flow 
regions in this study. 

The boundary-layer equations shown above were solved using the relaxation 
method given by Van Dalsem and Steger (1985). The numerical algorithm is writ- 
ten to solve the equations in a transformed domain which is a uniformly spaced, 
rectangular coordinate region, using upwind differencing in the streamwise and 
spanwise directions, and central differencing in the normal direction. Details of the 
transformation and the numerical schemes used can be obtained from Van Dalsem 
and Steger (1985). 

In analyzing three-dimensional boundary-layer solutions, it is useful to define 
some dimensionless parameters such as the wall shear stress components t wx , T wy ; 
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and integral thicknesses 6 X , S y as 

du 

TyoX = IU 

oz 


dv 

t wxj = ZZZ u> 
OZ 



in which the subscript e refers to values at the edge of the boundary layer. Note 
that in these definitions, the normal coordinate z is scaled by Rez . 

Unlike two-dimensional flows, the integral thickness 6 X is not the correct mea- 
sure of the displacement effect of a three-dimensional boundary layer. For three- 
dimensional flows, Moore (1953) showed that the displacement thickness 6* has to 
account for cross-flow effects. He derived a partial differential equation for the dis- 
placement thickness, which, in terms of the integral thicknesses defined in equation 
(4.9), is given as 

d d 

ju,((S — b x )j ” ,. (ve^ — Ugbx/) — 0 (4.10) 

ox ay 

From a known solution of the boundary-layer variables, it is relatively easy to 
numerically solve the partial differential equation given by equation (4.10) and 
obtain the displacement thickness. For the geometries considered in this study, the 
flow was assumed to be two-dimensional at the inflow boundary, and hence, 6* = 6 X 
there. With known values of u e , v e , 6 X , and 6 y , an integration algorithm similar to 
that used for integrating the vorticity in § 11.E.3 was used to obtain the displacement 


thickness. 
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2. Two-dimensional du a l potential equations 

It is, of course, possible to use the three-dimensional formulation to compute 
two-dimensional flows. However, it is more convenient to present the dual potential 
equations for the two-dimensional case separately. 


a. Governing equations To retain the terminology used earlier, the normal 
coordinate is denoted as z instead of the traditional y. With this convention, 
the governing Navier-Stokes equations for steady, incompressible, laminar two- 
dimensional flow are: 

Continuity: 

u x + w z = 0 (4.11) 


x- momentum: 


p(uu x f wu z ) 


' ke (Ul 


+ Uzz) 


(4.12) 


e-momentum: 

p{uw x + ww z ) = -pz + ^(wn + w 2z ) (4.13) 

Using the dual potential decomposition defined in terms of the scalar potential 
function d> and the stream function \ as 


u — x Xz 


w = 4>z + Xi 


(4.14) 


the continuity equation is written as 


V 2 4> = o 

and with the vorticity defined as 


(4.15) 


{jJ — u z - w x 


(4.16) 
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the stream function equation is 

V 2 X = -w (4.17) 

The momentum equations are written in terms of the vorticity as 

p(uou' x + wuj z ) ~ J^i^xx + w 2r ) (4.18) 

Note the absence of the vorticity stretching terms in the two-dimensional vorticity 
transport equation. 

b. Boundary conditions Boundary conditions at any boundary surface are 
no-slip and tangency, that is, 

u = v = 0 (4-19) 

No-slip is enforced implicitly by evaluating the vorticity at the wall appropri- 
ately. Tangency is enforced by selecting the stream function to be a constant along 
the surface and computing appropriately. 

= constant; n • V 4> — 0 (4.20) 

As in the three-dimensional case, vorticity for the interaction scheme was ob- 
tained from a boundary-layer solution. The algorithm used to solve the boundary- 
layer equations was a two-dimensional version of that used for the three-dimensional 
case (Van Dalsem and Steger, 1985) and is described in Appendix B. 

B. Dual Potential - Boundary Layer Interaction Scheme 

As mentioned earlier, the present interaction scheme involves obtaining the vor- 
ticity from the boundary-layer equations. To solve the boundary-layer equations, it 
is necessary to determine the pressure from the dual potential solution. A schematic 
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view of this coupling method is given in Figure 17. Such an interaction scheme was 
proposed and implemented for two-dimensional flows by Steger and Van Dalsem 
(1985). 

1. Sol utio n strat egy for the int eraction s che me 

An iterative solution strategy, similar to the one used for inviscid flows is used 
for the dual potential - boundary-layer interaction scheme and is described below. 

At each grid point an initial guess is made for values of \ . v> and for 
and 1U3. The following iteration scheme is then implemented: 

]) Vect or poten tial fu n cti ons, ^,\,t p. Individually update th \ , and ij' from 
the Poisson equations given by equation (2.45) for assumed values of u>i,u;2, and 
u>3. Boundary conditions for these equations are again kept compatible with the 
consistency relation, equation (2.50). 

2) Scalar potential functi o n, <£ . Update <f> from the continuity equation (2.4 4) 
using the ADI-like algorithm for the transonic full potential function and previously 
updated values of i?,x, and xjj. 

3) Pressure, p . With the velocity components known, pressure is calculated 
by integrating the normal momentum equation from the far field (uniform flow, 
p = poo) to the wall. The numerical scheme employed is similar to the one described 
in §11. E. 3 for integrating the vorticity consistency equation. 

4) Vorticities, uq.u;?, and u>a. With pressure known, the vorticity components 
are evaluated from the solution of the three-dimensional boundary-layer equations 
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i 

u>i = -Vz Rei 

u>2 = Uz Re 2 (4.21) 

^3 ~ ~ My 

Nole that the bars on the variables at the right hand side of equation (4.21) refer 
to boundary-layer quantities. Also, in accordance with the boundary-layer assump- 
tion, the contribution of the normal velocity is neglected in evaluating wj and u>2- 
Furthermore, the scale factor Re 5 appears on the right-hand-side of equation (4.21) 
since the normal coordinate z is scaled by this quantity in the boundary-layer equa- 
tions. 

5) Test for convergence and, if necessary, return to (l). 

The numerical algorithm used to integrate the normal momentum equation is 
described below. 

2. Int egration scheme for pressure 

In the dual potential - boundary layer interaction procedure described above, 
the boundary-layer equations require the knowledge of pressure from the dual po- 
tential solution. However, pressure is not a dependent variable in the dual potential 
formulation and a method for calculating the pressure from the dual potential solu- 
tion needs to be devised. If the flow above the boundary-layer is irrotational, it is, 
in principle, possible to determine the edge of the boundary-layer (defined as the 
first location in the normal direction z where ^ > c; e = 0.9995 in this study) and 
evaluate the pressure from Bernoulli’s equation applied to the dual potential solu- 
tion at this height. However, this procedure was found to be inaccurate for various 
reasons. In particular, in regions of separated flow, there is considerable growth of 
the boundary-layer and the normal pressure gradient in this region is not necessarily 


small in the dual potential solution and it was found necessary to use an alternative 
method to determine the pressure. Furthermore, with the stretched grids used in 
the calculations, some form of interpolation may be required in determining the 
edge of the boundary layer. Using the normal momentum equation appears to give 
a reasonable estimate of the pressure near the wall from the dual potential solution 
for the geometries and flow conditions considered in this study. With the use of 
the normal momentum equation, interpolation is avoided, and it is also possible to 
compute the pressure for the case in which the outer flow is rotational. 

With 2 assumed to be in roughly the normal direction to the surface and ne- 
glecting viscous terms, the 2-momentum equation (2.2) is written as 

p{uw x -I- vwy -I- ww z ) = ~p z (4.22) 

In transformed coordinates, this equation becomes 

tzPc + PzPf, + UP< = -p{Uwc + Vwrj + Wu^) (4.23) 

This equation has the same form as equation (2.71) where / = p, a\ — £ z , a<i = 
r \ z , 0.3 — and 04 is the right hand side of equation (4.23). With the pressure 
known in the far field, that is, at / = LM AX, p = poo, equation (4.23) can be 
integrated in the — f-direction to obtain the wall pressure using algorithms similar 
to that given by equations (2.73) - (2.74). First-order backward differencing is used 
at the first f = constant plane, and second-order differencing at all subsequent 
planes. 

jb Interaction scheme for separated flow 

For attached flows, the solution strategy given above is easy and straightforward 
to implement. For separated flow, it is not possible to use pressure as the forcing 
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function, and a way of calculating or updating the shear stress needs to be devised. 
Van Dalsem and Steger (1985) found the following scheme to be fast and reliable 
for updating the shear stress: 

* 1 =Cx<“'(p”-lC l ) (<•«) 


and used a value of the relaxation parameter w % 10. In equation (4.24), p w refers 
to the pressure at the wall which is obtained from the dual potential solution by 
integrating the normal momentum equation. Thus, they updated the shear stress 
from the difference between the viscous and the inviscid pressures at the wall. The 
viscous pressure p was obtained by using Bernoulli’s equation at the edge of the 
boundary layer as 

p = ^(i - («t 4 *’?)! ( 4 - 25 ) 


From equation (4.8), it can be seen that there is a direct relation between the 
shear stress and the pressure gradient. In an iterative solution method, this relation 
can be used to devise a new scheme to update the shear stress. Equation (4.8) can 
be rewritten as 


r = 






U 


n+1 


u 


rH 1 


+ T. 


n+1 


^2 


(4.26) 


When the interaction scheme between the dual potential and the boundary-layer 
formulation converges, the residual term r n on the right-hand-side of equation (4.26) 
should be zero. Hence, a formula for updating the shear stress can be formed as 


T 


n+1 

wx 


— uir 


(4.27) 


where u is a relaxation parameter. Note that the scheme for updating the shear 
stress at the wall given by equation (4.26) - (4.27) does not require the assumption 
of irrotationality at the edge of the boundary-layer. In fact, it does not require 
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any assumptions or values of variables at the edge of the boundary-layer provided 
the pressure gradient from the dual potential formulation can be obtained without 
the knowledge of any edge values. This is possible when the normal momentum 
equation is used to obtain the pressure from the dual potential solution. 

C. Flow Over a Two-Dimensional Trough 

The dual potential - boundary layer interaction procedure described above was 
used to compute the flow over a two-dimensional trough which was first studied 
by Carter and Wornum (1975), and subsequently by Kwon and Pletcher (1981), 
Veldman (1981) and Edwards and Carter (1985). The algorithm used to solve the 
boundary-layer equations was a two-dimensional version of the three-dimensional 
one given by Van Dalsem and Steger (1985) and is given in Appendix B. The bound- 
ary layer algorithm has been validated for a number of two- and three-dimensional 
flows by Van Dalsem and Steger (1985, 1986a). In the present study, the boundary- 
layer code was again verified by computing flow over a flat plate, both in the direct 
and inverse modes. In all cases, the agreement with known analytical /numerical 
results was good. The dual potential procedure was verified in Chapter III for in- 
viscid flows. Details of the numerical algorithms for the two-dimensional case are 
given in Appendix B. The dual potential - boundary layer interaction procedure 
was verified by computing the flow over a flat plate and a detailed description of 
these results will not be presented here. 

1. Geometry and grid 

The geometry considered here is shown in Figure 18. It consists of a flat plate 
with a trough located at some distance from the sharp leading edge. The surface 
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of this geometry is described by the function 

s(x) = -f sech(4x - 10) (4.28) 

The reference length for this configuration is L — 1 m. The free stream Reynolds 
number based on the reference length is 80,000. In the present calculations, the 
inflow and outflow boundaries for the dual potential set of equations were set to x = 
— 2.5 and x = 7.5, respectively and the outer boundary was set at z - 5 with a mesh 
of 219 x 81 grid points in the x and z directions, respectively. Clustering functions 
given by Vinokur (1983) were used in order to cluster points in the streamwise 
direction at the location corresponding to the maximum depth of the trough, that 
is, x — 2.5. A stretching function was used in the normal direction to obtain a 
fine grid in the vicinity of the wall, with the minimum spacing at the wall set to 
0.0001 . The grid used to obtain the dual potential solution is shown in Figure 
19. The interaction region extended from x — 1.0 to x = 4.0 which correspond to 
121 points in the streamwise direction. The boundary-layer grid extended up to a 
distance of z = 0.2 with 55 points in the normal direction. The same distribution 
of points in the normal direction as in the dual potential grid was used, hence the 
need for interpolating the vorticity was avoided. For the boundary-layer, a Falkner- 
Skan solution for zero pressure gradient scaled to the appropriate location in the 
x-direction (in this case, x — 1.0 from the leading edge of a flat plate) was used as 
the inflow boundary condition. 

Once again, body fitted curvilinear coordinates were employed, and the flow 
domain mapped to a uniformly spaced rectangular coordinate region. The details 
of these transformations are very similar to that given for the three-dimensional 
transformation used in Chapter II and are given in Appendix B. 


80 



Figure 19. View of grid for the two-dimensional trough, t — 0.03 
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2. Results and discussion 

To verify the boundary-layer code for separated flow, solutions were obtained 
from the code separately with the inverse forcing function specified. The shear stress 
distribution used as the inverse forcing function was obtained from the results of 
Davis et al. (1986) who used the vorticity /stream-function formulation to solve the 
Navier-Stokes equations. Figure 20 is a plot of the wall shear stress r wx which was 
prescribed for this case. Figure 21 show comparisons of the pressure p — 5(1 - uj) 
where u e is the velocity at the edge of the boundary layer, with the results of 
Veldman (1981) who studied the same configuration using a quasi-simultaneous 
viscous-inviscid interaction procedure. The present results are compared with those 
of Veldman because Davis et al. (1986) did not present values of the pressure in 
their report. The figure shows good agreement between the present boundary-layer 
results and that of Veldman. Convergence was obtained in about 100 iterations 
requiring 20 seconds on a CRAY XMP computer. 

Attempts to compute the same flow with the dual potential - boundary layer 
interaction procedure failed, however, as it was not possible to obtain a converged 
solution. 

In order to examine the reason for this failure, the vorticity was computed from 
the converged boundary-layer solution (for a specified wall shear stress distribution) 
and used to obtain a converged dual potential solution. Results for this case are 
presented in Figure 22 which is a plot of the streamwise variation of the pressure. 
The pressure from the dual potential solution was calculated by integrating the z- 
momentum equation from the far-field down to the vicinity of the wall as described 
earlier. Figure 22 is a comparison of the pressure between the dual potential so- 
lution and the result of Veldman (1981) for this case. It can be seen that there is 




Figure 21. Pressure distribution for the two-dimensional trough from the 
boundary-layer solution for a specified wall shear stress distri- 
bution, / — 0.03 
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Figure 22. 


Pressure distribution for the two-dimensional trough from the 
dual potential solution for specified vorticity field 
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considerable disagreement between the two curves in the vicinity of the separated 
flow region. This discrepancy between the two solutions is the apparent reason for 
being unable to obtain a converged interaction solution for separated flows. 

In general, the boundary-layer equations are sensitive to the forcing functions, 
that is, pressure or w'all shear stress (cf. Klineberg and Steger, 1974; Lee and 
Pletcher, 1986). Any error in these functions generates an erroneous vorticity solu- 
tion which in turn feeds back into the dual potential solution, in this case causing 
an unstable interaction process. 

It should be pointed out that the interaction procedure described in this study 
is different from the more conventional viscous-inviscid interaction procedures in 
which boundary layer effects are imposed on the outer flow as a correction to the 
shape of the body surface. In the current procedure, the dual potential solutions 
are solved over the entire flow domain including the region where viscous effects 
are important, and the viscous effects are accounted for by injecting vorticity from 
the boundary layer equations. Since, in principle, the only approximation made 
is in using the boundary-layer equations rather than equation (4.18) to provide 
the vorticity, the issue arises as to whether this vorticity is an accurate input for 
the governing dual potential equations. To answer this question, the dual poten- 
tial formulation was modified by obtaining the vorticity directly from the vorticity 
transport equation (4.18) thus avoiding any boundary-layer assumption. To be 
consistent with previous data, however, the specified shear stress described earlier 
was used to prescribe the wall vorticity. Thus, even though the full Navier-Stokes 
equations are solved for this case, the calculation used as its vorticity boundary 
condition the previous known solution obtained by other means. The numerical 
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algorithm for solving equation (4.18) is similar to the one described in §11. E. 2 and 
is given in Appendix B. 

In practice, the boundary-layer vorticity was used at the first two points in the 
normal direction, that is, at / = 1,2, to avoid any of the numerical errors at the 
walls due to the stretched grids (cf. Appendix A). Solutions obtained from the 
Navier-Stokes formulation for this case are presented in Figure 23, which is a plot 
of the pressure from the dual potential solution along the length of the trough, 
compared with the result obtained by Veldman (1981). The agreement between the 
two pressures is noticeably better from the case where the boundary-layer vorticity 
was used in the viscous layer (Figure 22). To further study this case, vorticity 
profiles from the bound ary- layer and Navier-Stokes solutions are compared at select 
locations in the streamwise direction in Figures 24-26. From these figures, it can 
be seen that there is considerable discrepancy between the two vorticities in the 
separated flow region. 

It should be pointed out that at the flow conditions for this case, namely the 
moderate Reynolds number and the depth of the trough, the bound ary- layer is 
relatively thick (almost twice the depth of the trough) and there is considerable 
variation of the pressure gradient in the normal direction. This fact is illustrated 
in Figure 27 by plotting pressure from the inviscid solution at the wall and at a 
location corresponding to the edge of the boundary layer for this configuration. 
Figure 27 shows that there is significant variation of the pressure gradient in the 
normal direction within the shear region. Hence, the boundary-layer assumption is 
questionable in such regions and the type of interaction considered here is perhaps 
not appropriate for this case. A more suitable form of interaction is needed in which 
the boundary-layer equations are solved in regions very close to the wall (e.g., for the 
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Figure 23. Pressure distribution from the Navier-Stokes solution with vor- 
ticity specified at the first two points from the wall, 1 = 2 
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Figure 24. Comparison of vorticity profiles between Navier-Stokes 

boundary-layer solutions, Navier-Stokes; 

boundary-layer 
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Figure 25. Comparison of vorticity profiles between Navier-Stokes and 
boundary-layer solutions, a) x = 1.62; b) x = 2.41 
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first twenty points considered in this study), with appropriate vorticity boundary 
conditions derived from the dual potential solution applied at the “edge” of the 
boundary-layer grid (instead of the usual irrotational assumption). Such a scheme 
could perhaps form the basis of a detailed study on the breakdown of the boundary- 
layer equations in separated flow regions. 

From the results presented in this section, it is possible to conclude that the 
dual potential formulation of the Navier-Stokes equations can be applied to calcu- 
late separated flows, but caution has to be exercised in using the boundary-layer 
equations to supply the vorticity in regions away from the wall. 

D. Flow Over a Three-Dimensional Trough 

The dual potential - boundary layer interaction procedure was applied to com- 
pute the incompressible laminar attached flow over a three-dimensional trough con- 
figuration and results for this case are presented in this section. 

L Grid and geometry 

The geometry considered here is a three-dimensional version of the geome- 
try considered in the previous section and a schematic view is given in Figure 
28. Edwards (1986) obtained a solution for this geometry using an interacting 
boundary-layer algorithm coupled with an inviscid small-disturbance analysis, while 
Davis et al. (1986) obtained a viscous flow solution for this geometry using a 
vorticity /stream-function type formulation for the Navier-Stokes equations. The 
surface of this geometry is described by the function 


z{x,y) = — tsech(4x - 10)sech(4y — 6) 


(4.29) 
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Once again, the reference length for this geometry is L = 1 m. The shape of the 
trough is symmetrical about the y — 1.5 plane. Results for inviscid, and viscous 
flow at a free stream Reynolds number based on the reference length of 8000 were 
compared with those of Davis et al. (1986) and Edwards (1986). These results 
show that the present dual potential procedure is capable of accurately predicting 
three-dimensional inviscid and viscous laminar flows. 

Figure 29 shows a view of the computational grid in the x-z plane which was used 
to obtain a three-dimensional dual potential solution. The x- and ^-computational 
boundaries were located in the region 1 < x < 4, 0 < y < 3. A total of 41 x 40 / 51 
grid points were used. A uniform grid distribution was used in the x- and y- 
directions. The 51 points in the 2 -direction were distributed from the lower surface 
to 2 = 5 with a minimum spacing at wall equal to 0.0001. In all calculations, 
periodic boundary conditions were used in the y-direction. For inviscid flow, the 
following boundary conditions were used at the other boundaries: 

Inflow : 

x = x/, = d m 0; <t> x = Uoo (4.30) 

Outflow : 

XxX = V'zx = ^ XX = 0; 4>x — U OC (4.31) 

Wall : 

0 = X = 0> = -(tfz + Xy), <t>n = o (4.32) 

Far field : 

= X = 0, V>z = ~(0x + Xy), <t> = UooX (4.33) 

Note that the flow at the inflow boundary (i = 1) is assumed to be uniform because 
the trough thickness asymptotically goes to zero away from the center so that at this 
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Figure 29. View of the sheared grid for the three-dimensional trough in the 
x-z plane at y — 1.5 
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boundary, the thickness is practically negligible ( t = 0.00016). Another reason for 
choosing this boundary condition was to match those used by Davis et, al. (1986). 
It will be shown that for the conditions studied, these boundaries are adequate to 
capture any elliptic effects, and provide a smooth solution for the flow field. 

For viscous flow, the same boundary conditions were used at all but the inflow 
boundary, where \ was computed by solving a one-dimensional Poisson equation 
for a prescribed vorticity. The vorticity for this boundary was obtained from the 
boundary-layer profile for the velocity at this location, which was obtained from a 
Falkner-Skan solution (e.g., White, 1974). 

Figure 30 shows the grid in the plane of symmetry which was used to obtain the 
boundary-layer solution. The edge of the grid in the ^-direction is not uniform in 
x but varies as the square root of x. At the inflow boundary, the grid extends from 
z — 0 to z — Zmax = where 6* is the displacement thickness for a flat plate 
boundary-layer at a distance of x 1 from the leading edge, and o is a coefficient to 
account for boundary-layer growth due to adverse pressure gradients. In this study, 
a value of a = 50 was used giving a value of z m ax — 0.2 at the inflow boundary. 
Once again, a nonuniform grid was used in the z-direction, with a total of 51 points 
extending from z — 0 to z = z ma x with a minimum spacing at the wall of 0.0001. 

Since the grids in the z-direction in the boundary-layer and dual potential codes 
are not necessarily identical, a cubic spline function interpolation scheme is used 
to interpolate the vorticity components from the boundary-layer grid to the dual 
potential grid. For the geometries considered in this study, the variations in the x- 
and ^-directions between the two grids is small and vorticity was interpolated only 
as a function of the normal direction z (or z). 
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2. Inviscid results 

As another simple verification of the numerical procedure, a three-dimensional 
inviscid solution was calculated for a f = 0.06 trough with the scalar potential 
function alone. Figures 31 and 32 are contour plots of the u- and u-velocity distri- 
butions on the surface of the trough. These figures show that the predicted velocity 
distributions are smooth and that the flow becomes almost uniform at the outflow 
boundaries. The three-dimensional flow produced by the trough is illustrated by 
Figure 32, which shows that upstream of the maximum depth of the trough, the 
pressure gradients created by the increasing depth of the trough push the flow to- 
wards the plane of symmetry, while in the downstream half of the geometry, the 
flow is in the direction away from the plane of symmetry. More detailed plots of 
the x- and y- components of the velocity are given in Figures 33 and 34 at two 
different y-planes and compared with ot her numerical results. Figure 33 is a plot of 
the streamwise velocity u on the surface of the trough versus x , while Figure 34 is a 
plot of the spanwise velocity v on the surface of the trough versus x. Comparisons 
with other numerical computations as taken from Davis et al. (1986) are shown in 
these two figures. The agreement with the results of Edwards (1986) and Davis et 
al. (1986) is relatively good for the u velocity with the present results having values 
in between these other results. However, for the v velocity (Figure 34), the present 
results agree well with those of Edwards and differ slightly from those of Davis et 
al. 

It should be pointed out that the results of Davis et al. (1986) were obtained on 
a81 x41 x81 grid. In the present calculations, which were performed on a CRAY 
XMP-48 with a maximum allowable memory of 4 million bytes, it was not possible 
to obtain solutions with finer grids. Hence, results were obtained with coarser grids 
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Figure 32. Contour plots of inviscid v velocity at the surface, t — 0.06 
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to determine the effect of grid size on the solutions. Figures 35 and 36 show results 
obtained with a 41 x 20 x 51 grid compared with those from a 41 x 40 x 51 grid. In 
these figures, the surface velocity components are plotted versus x at two spanwise 
locations for the two different grids. The figures show little or no difference in the 
results due to halving the number of grid points in the y-direction. All subsequent 
calculations, unless otherwise indicated, were performed on a 41 x 40 x 51 grid 

The present calculations required approximately 900 iterations and took about 
300 seconds on the CRAY-XMP computer. 

3. Viscous results 

A three-dimensional viscous solution was calculated for a t = 0.03 trough at 
an upstream reference Reynolds number of 8000 with the interaction procedure 
described in §IY.B. For this configuration, the flow remained attached over the 
entire surface. The flow was once again assumed to be two-dimensional at the 
upstream boundary. The boundary layer velocity profile was found from a solution 
of the Falkner-Skan equations for zero pressure gradient flow (e.g., White, 1974) 
and scaled to the appropriate location in the x-direction (in this case, x = 1.0 from 
the leading edge of a flat plate). 

As a measure of the agreement between the dual potential and boundary-layer 
solutions, the pressures from both schemes (that is, the pressure which is specified 
to the boundary-layer equations, and the pressure from the boundary-layer solution 
which is calculated from Bernoullis’s equation applied at the edge of the boundary- 
layer) at the plane of symmetry are presented in Figure 37. The agreement between 
the two is relatively good indicating the convergence of the interaction scheme. Some 
discrepancy is attributed to the different formulations and numerical algorithms. 
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Figures 38-43 are contour plots of the wall shear stresses, the integral thick- 
nesses, the displacement thickness, and the boundary-layer pressure, respectively, 
for the three-dimensional trough geometry. These figures give a qualitative picture 
of the flowfield. The growth of the boundary layer is evident from the displacement 
thickness plots. It can be seen from these figures that the flow starts as two- 
dimensional, and then following the geometry, becomes three-dimensional towards 
the center of the trough. However, the flow does not become two-dimensional as it 
approaches the outflow boundary, even though the geometry reverts to a flat sur- 
face, but reaches a quasi three-dimensional state. This is clearly evident in Figure 
42, which is a contour plot of the displacement thickness. A comparison of Figures 
40 and 42 also illustrates the effect of the cross-flow on the displacement thickness 
for the trough geometry. The presence of cross-flow has decreased the displacement 
effect of the boundary-layer from what it would have been if the flow was purely 
two-dimensional. 

A qualitative estimate of the viscous effects on the inviscid flow is obtained by 
comparing the pressure from the interaction solution in Figure 43 with a contour 
plot of the inviscid pressure for this case which is given in Figure 44. The two 
figures show that even for this attached flow case, there is considerable interaction 
of the boundary-layer with the outer inviscid flow. In particular, the pressure from 
the interaction solution varies from -0.010 to 0.025 whereas the inviscid pressure 
varies from -0.02 to 0.08 in the flow domain. 

Detailed results for the boundary-layer parameters are presented in the line 
plots given in Figures 45-49. These figures are plots of the streamwise variation of 
the various flow parameters at different span locations. Since the flow is symmetric 
about the y = 1.5 plane, results on only one side of the plane of symmetry are 
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Figure 38. Contour plot of x-component of wall shear stress, t = 0.03 



Figure 39. Contour plot of y-component of wall shear stress, t = 0.03 



Figure 40. Contour plot of integral thickness in the i-direction, t = 0.03 






Figure 41. Contour plot of integral thickness in the y-direction, t = 0.03 




Figure 42. Contour plot of displacement thickness, t - 0.03 
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Figure 43. Contour plot of pressure from the interaction solution, t = 
0.03 
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plotted. Comparison with other numerical results as obtained from Davis et al. 
(1986) are presented in Figures 50 and 51. Figure 50 is a plot of the variation 
of t wx in the streamwise direction at two different y locations. The agreement 
with other results is relatively good. Figure 51 is a plot of r wy versus x again at 
two different y locations. In this figure, the agreement with the small-disturbance 
interaction results of Edwards (1986) is very good, while those of Davis et al. (1986) 
differ considerably. Davis et al. (1986) believed that the discrepancy in the two 
results was caused by the grid spacing, though no conclusive studies are available 
for explaining the discrepancy. It should be noted that it is difficult to study the 
effect of grid spacing in three-dimensional problems even with the supercomputers 
available today. As for the inviscid case, grid effects were examined by halving the 
number of points in the y direction. Once again, there was little or no difference in 
the results, and hence detailed comparisons are not presented here. 

Further comparisons between the two solutions (boundary-layer and dual po- 
tential) are presented by plotting the velocity profiles from both the dual potential 
and the boundary layer solutions in Figures 52 through 54. Figure 52 is a plot of the 
u velocity profiles in the plane of symmetry at different x locations. The growth of 
the boundary layer and the effect of the adverse pressure gradient near the center of 
the trough are visible in this figure. The dual potential and boundary-layer profiles 
show very good agreement. Detailed plots are given in Figures 53 and 54 which are 
plots of the velocity profiles in the maximum trough depth region at two different 
y planes. Again the agreement between the two solutions is relatively good. The 
no-slip boundary condition is not imposed on the vector potential functions directly 
and there is a small slip velocity in the dual potential profiles of the order of 2% of 
the free stream velocity. This error is partially geometric in origin, produced by the 
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Figure 50. Comparison of the x-component of wall shear stress, t = 0.03. 
a) y = 1.5; b) y = 1.2 
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Figure 52. Velocity profiles from the dual potential and boundary-layer so- 
lutions at plane of symmetry (y = 1.5), t = 0.03, dual 

potential; boundary-layer 
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Figure 53. Comparison of u-velocity profiles from the dual potential and 
boundary-layer solutions at y = 1.5, t = 0.03, dual poten- 
tial: boundary-layer 
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Figure 54. Comparison of u-velocity profiles from the dual potential and 
boundary-layer solutions at y — 0.75, t = 0.03, dual po- 
tential; boundary-layer 
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highly stretched grids which need to be used in the 2-direction. A one-dimensional 
analysis of the cause of this error is given in Appendix A. Another possible source 
of error is interpolation errors due to the different grids used in the two schemes. 

The three-dimensional interaction calculations required approximately 1400 it- 
erations which took about 2500 seconds on the CRAY XMP. The iterative procedure 
was terminated when the residuals for both the dual potential and the boundary- 
layer algorithms had dropped by at least six orders of magnitude. As a practical 
convergence criteria of this method, the calculation was assumed to be converged 
when the relative changes in the pressures and the wall shear stress was less than 
10~ 4 between successive iterations. The calculations were performed in the direct 
mode for the boundary-layer algorithm. Essentially the same solutions were ob- 
tained by using the inverse mode. However, the calculations in the inverse mode 
required more iterations. 

The results described in this section show that the dual potential formulation, 
when coupled with an appropriate boundary-layer algorithm, can be used for solving 
three-dimensional viscous flow problems in which the viscous layer remains thin as 
in the case of attached flow. 



V. CONCLUSIONS 


A. Concluding Remarks 

A dual potential decomposition scheme has been devised for computing three- 
dimensional rotational flows. With this formulation, rotational effects can be stud- 
ied by adding the appropriate vorticity into the governing equations for the vector 
potential functions. 

In the first part of this study, a finite-difference procedure for solving the dual 
potential equations has been developed and applied for calculating the inviscid 
but rotational flow through indraft wind tunnels, including the effects of screens 
and vanes. Actuator disk theory was used to model the stagnation pressure loss 
produced by flow straightening devices such as screens and vanes which are located 
at the inlet entrance. In particular, the flow through the 80- by 120-ft wind tunnel 
at the NASA Ames Research Center was successfully simulated using semi-empirical 
loss-coefficients to account for the pressure drop across the screens and vanes and 
the turning effect of these devices. The procedure was applied to compute the flow 
for two different inlet vane and screen configurations. The numerical solutions were 
in satisfactory agreement with experimental and other numerical results. With 
the use of this procedure, the relative effect of various passage geometries, flow 
straightening devices, and external wind conditions on test-section flow quality can 
be studied. 

In the second part, viscous effects were modeled by adding vorticity in the 
boundary-layer region. This vorticity was obtained through a suitable boundary- 
layer interaction procedure. The numerical results obtained for attached flow over 
a three-dimensional trough compared favorably with other numerical solutions. For 
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separated flow, the numerical results indicate that the present formulation of the 
dual potential - boundary layer interaction procedure fails to converge. For two- 
dimensional separated flows, the dual potential form of the Navier-Stokes equations 
was solved using a prescribed vorticity boundary condition. A comparison of the 
vorticities from the boundary-layer and the Navier-Stokes solutions indicates that 
for the type of interaction considered in this study, the boundary-layer solution does 
not. give a representation of the vorticity which is consistent with the dual potential 
formulation in regions where the boundary-layer is relatively thick as in separated 
flow. 

Overall, a versatile and flexible tool has been developed by which it is possible to 
simulate a wide variety of flows with little or no changes to the numerical algorithm. 
In particular, it is possible to study rotational effects by adding the appropriate 
vorticity into the governing dual potential formulation. 

B. Recommendations for Future Study 

The dual potential velocity decomposition has been shown to be a viable alter- 
native to primitive variable formulations for solving many three-dimensional flow 
problems. In this study, the formulation was used to develop a procedure for com- 
puting both inviscid and viscous flows with relatively small changes in the algo- 
rithms. 

The dual potential formulation has already been applied successfully, in two di- 
mensions to compute the inviscid transonic flow over airfoils (Chaderjian and Steger, 
1985). The extension to transonic flow in three dimensions seems straightforward 
and needs to be implemented. While the present procedure has been successfully 
used to compute three-dimensional attached viscous flow, and some of the problems 



associated with separated flow have been addressed in this study, the issue of sepa- 
rated flow remains an unsolved problem. In the interaction procedure used in this 
study, the boundary-layer equations were solved in the primitive variable form. It 
would perhaps be advantageous to use a vorticity /stream-function formulation of 
the boundary-layer equations in an alternate interaction procedure. 

The dual potential formulation has not been clearly understood or developed for 
unsteady, compressible flows, even in two dimensions. Hence, further w'ork remains 
to be done in extending the formulation to these more complex flows. 

Finally, the treatment of turbulent flows is still one of the unsolved problems in 
fluid mechanics today. The use of the dual potential formulation for investigating 
turbulent flows needs to be examined; in particular, the selection of vorticity as the 
dependent variable in this formulation may have some advantages in this regard. 
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VIII. APPENDIX A: GEOMETRIC ERRORS 


In this appendix, the geometric errors introduced by using stretched grids in 
the computational domain will be analyzed by considering a one-dimensional model 
problem of the vorticity /stream-function case. 

Consider Laplace’s equation in one dimension 

Ipyy = Oj U — Ipy (A.l) 

solved on a simple geometrically stretched grid defined by 

dy{k) = ( k ~ l dyo; k = 1, . . . , KM AX - 1 (A.2) 

where dy(k) = y(k + l) — y{k). and KM AX is the total number of points, xj) = y 
is a solution of this equation for the boundary conditions 

at y = 0 : %l> = 0; at y = ymax : xjjy = 0 (A.3) 


When the domain is transformed into a uniformly spaced computational domain 
defined by y — y(y ), the governing equation (A.l) is transformed as 


dy 


\J dy J 


= 0 


(A.4) 


where J is the Jacobian of the transformation defined as 


J — 1/j/n 


(A.5) 


and 


r]y = 1/ifcj 

Hjscaa no 


(A.6) 
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In the discretized domain, the left-hand-side (LHS) of equation (A. 4) is approxi- 
mated as 


LHS = \{ yr , k + y^ +1 )(4- + - V^) 


Vr i k Vr >k + 1 




4 w ''fc ""'fc- 1 ' ' u i ,/2 

% y« k _ l 


){^k - V’fc-i) 


(A.7) 


The metric quantities for this case can be calculated analytically for the particular 
stretching form that is used as shown below 


k y 

1 yi = 0 

2 V2 = dyo 

3 y 3 = dyo( 1 + t) 

4 y 4 = dyo{ 1 ~ < 4 c 2 ) 

5 y 3 = dyo( 1 + f -I- £ 2 + e 3 ) 


yr, 

yrn = y2 - yi = rf y° 
y^ = 2(^3 - yi) = jcM 1 + f ) 

yr? 3 = \(V4- V2) = \td.yo{ 1 + <) 
Vt]4 = \{V5 ■- ys) = 2e 2 dyo(l + t) 


(A.8) 


Setting 4' — y an d substituting the various terms in equation (A.7) for different 
values of k the left-hand-side can be evaluated exactly. For k — 3,4,..., KM AX - 2 
it can be shown that LHS = 0, and hence the analytical solution satisfies the numer- 
ical approximation exactly for any values of the stretching parameter e. However, 
for k - 2, and for k = KM AX — 1, LHS is not zero but is evaluated as 


lhs 1=2 = 


1 

— £ 

2 



1 (3 + c)(f 2 + 2c + 5) 
8 (1 + c 2 ) 


(A.9) 


In this case, LHS — 0 if, and only if, c = 1 which corresponds to a uniform mesh. 

The geometric errors introduced in the computation domain for the case of the 
simple geometric stretching have been analyzed above. The exact nature of the 
error would depend on the type of grid used. This kind of error has been analyzed 
by others ( cf . Hindman, 1981; Thomas and Lombard, 1978). It has been shown that 
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the error is localized to the boundary points for the kind of stretching considered 
here. This error is demonstrated in a solution of the one-dimensional homogeneous 
problem which is presented in Figure Al. Figure A1 is a plot of ‘velocity’ ( u = xjj y ) 
versus y for the case in which ymax — 5.0, KM AX - 81, dyo = 0.0001. From 
the figure it can be seen that the error is localized and does not affect the solution 
for the interior points. In particular, when the nonhomogeneous Poisson equation 
is solved for boundary-layer type vorticity, the error is manifested as a small slip 
velocity as shown in Figures A2 and A3. In this case, the equation tp yy = —u> is 
solved where u> is computed from a polynomial equation for Blasius flow with a 
boundary-layer thickness of 0.02. The slip velocity at the wall is of the order of 
1.5% for this case. 

One of the commonly used cures for this problem is the free stream subtrac- 
tion concept introduced by Pulliam and Steger (1980) whereby an exact analytical 
solution of the governing differential equation for uniform flow is subtracted from 
the right-hand-side of the numerical approximation. Unfortunately, this technique 
requires the knowledge of an exact analytical solution which may not always be 
known. 

Other solutions for this problem involve the computing of metrics in special 
ways or the use of special discretization formulas at boundaries (Vinokur, 1986). 
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IX. APPENDIX B: TWO-DIMENSIONAL EQUATIONS 

In this appendix, the governing equations, transformations and numerical algo- 
rithms for the two-dimensional incompressible dual potential and boundary-layer 
formulations which was used in §IV.C will be presented. The numerical algorithms 
are simplifications of the ones presented in Chapter III and hence some details may 
not be included here. 

A. Dual Potential Form of Equations 

1_. Governing equations 

For the dual potential decomposition given by 


u = 4> x -Xz-, w = 0* + x* (B.l) 

the governing equations are given as 

4>xx + <t>zz — 0 (B.2) 

X xx + Xzz = -w; w = u z - w x (B.3) 

uu>x T uju)z — — (wji "f ijJzz) (B.4) 


2. Transformation to general coordinates 

Using a general coordinate transformation defined by 

f = f(x,z), f = f(x,z) (B.5) 

in terms of the independent variables £, f , the Cartesian velocity components can be 
evaluated from the potential functions using the chain rule of partial differentiation 


UK 


on equation (B.l) as 

u = {Zx<t>S + $x<t>s) - UzX' + CzX?) 
w = + Cz<f>;) + Ux\c + ciXc) 

The metric quantities used in the above transformations are defined as 

Zx = JSs, 6 Jx c 

Cz = -Jzt, C 2 = Jic 

and J is the Jacobian of the transformation given by 


(B.6) 


(B-7) 




(B.8) 


Unsealed contravariant velocity components of the velocity vector q can be 
defined as 

U = i x u + ( z w 

(B.9) 

W = CzU + Cz«’ 


The contravariant velocity components are conveniently split between the scalar 
potential contribution and the vector potential contribution as 

U =r U 4, + U* 


w +W X 


(B.10) 


where 

U * = A“<t>c 4 
W * = A^<f>c + A* <f>, : 

U x = -Jx < 


(B.l 1 ) 


W x = J X C 

a^ = v^-ve = d + d 

= Vf ■ Vf = d + C* 

•A'* 5 ’ = • Vf = CxCx + £zCz 


(B.12) 
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3. Governing equations in transformed coordinates 


The continuity equation can now be written in transformed coordinates as 




+ 


{A^<f>c + A«4> ( ) 


= 0 


(B.13) 


The Poisson equations for the stream function can be transformed similarly em- 
ploying the same metric groupings and can be written as 


[A^X‘ + A&Xs) 


+ 


{A&X£ a A°'x, : ) 


while the vorticity transport equation (B.4) is written as 


(Uuc + Wu ( ) = 


Re 


' A^ujc + A&u s \ ( A^uc + A“ 


u; 

J 


u>- 


+ 


(B. 14) 


(B. 15) 


4. Numer ical algorith ms 

The continuity equation for 4>, equation (B.2), is updated using an approximate 
factorization (ADI-type) algorithm in delta-form. The differencing scheme in AF 
form is given by 

a£? aCC 

,n+l _ 


1 1 - hs^n ( i\] - 


= hu> 


_ u* - W* 1 

^('t r r + +6 ? (- T r-floo 


(B.16) 


The superscript n refers to the iteration level. The differencing of the contravariant 
velocity terms is illustrated by considering U ^ below 

U* +i = + \ {a%Mi+> + A ?Wj) ( B - 17 > 

The free-stream subtraction term f?oo appearing in the right-hand side of equa- 
tion (B.16), accounts for incomplete metric cancellation and is given by 




(B.18) 
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An ADI algorithm is used to solve equation (B.16). Rewriting equation (B.16) 
in the form 

LcL ? (4> n+1 - <j> n ) = R (B.19) 

the ADI algorithm is implemented in two steps as 

LcA4>* = R 

L ; A<£ n+1 = A<p* (B.20) 

<f> n+1 = 4> n + A 4> rt 

The algorithm given by equation (B.19) requires only a series of scalar, tridiagonal 
inversions and is, therefore, solved efficiently. 


The Poisson equations for the stream function is updated in a similar manner, 
which results in the AF form, 


C C 

i - hs,(^rk\\i - - x“) = + +s< <*r>" + j 


,A' : ' 


V 


w 


where 


X >+1 “ A T+l A tX] + o (‘ A 5+l i?X J+ 1 + A T 6 < X j) 


(B.21) 


(B.22) 


xY + l = A ]ii A <Xi + \ i + AfSnXi) 

The scheme given by equation (B.21) is implemented by the same two-step ADI 
algorithm given by equation (B.19). 

The differencing used with the vorticity transport equation is illustrated by 
considering equation (B.15): 


(Vu ( + Wu s ) = — 


+ A^Uc 


+ 


-f A^u< 


(B.23) 
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Assuming U is generally larger than W, three-point, second-order-accurate cen- 
tral differencing is used in f, and three-point, second-order-accurate upwind differ- 
encing is used in f. Upwind differencing of Use is automatically achieved using 


Uuc - + U~ b{<jj 


where 


(7+ = U + \ U \ 


and U — 


U - |17| 


2 2 

The convection equation is then solved using the AF relaxation algorithm 


(B.24) 


(B.25) 


/ + hU + 6 b c + hWS < - h^-L 

■» s Re s 


= -h + U-6* u n + W6(U n - ~ J 

where h is another relaxation parameter ( h > 0) and u/ and _<.’ W are defined as 

u/ i =~ A ss ! A cujj + - ^A" - • A' h- j.' ] 

•7+2 J+5 s 3 2 V J-i 1 ■ 1 ./ ■ J J 


(^r) ^ (' + «'-*/) K 1 -^) 


r+(.b..n , rr- cf, n , H/ r , n J It /. V \n . c / W > 


(B.26) 


= 


i - + 2 (-^?4 ] ^ W /-H + 


(B.27) 


Because central differencing is used in the ^-direction, numerical dissipation is added 
to the differencing scheme as 

7 / a C C \ ^ 1 

7 + + + hB'/> s . - fe— ^ - 3/i|Vy|AV|. : (/ + /tf/ 6 *) 1 - 

= -/> |t7 + ^w n + U~6[u> n + ~ [^(w r ) B + 6<(u W ) n ] + !W|(AV) 2 |, 

(B.28) 

At each £-plane, the algorithm requires a series of simple tridiagonal inversions in 
C- 


Tangency or no-flow-through is imposed on a boundary surface by setting the 
appropriate contravariant velocity component to zero. On a f = constant plane, 
W is set to zero. Tangency is enforced through implicit boundary conditions on 0, 
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which are obtained by solving the continuity equation at the half-cells neighboring a 
solid boundary. The center of such a cell is located at (y,/ + |). The finite-difference 
form of equation (B.13) can be written as 


- (l)]- 

A s c 


V' 




l+k 




Ac/ 2 


= 0 (B.29) 


Evaluating variables at / + \ by weighted averaging between / and / + 1 as. for 


example, 


1 


U J+Ll+^ - ~4 U J + y + + ~ (/ + 4 + 


(B.30) 


where 


V iHJ - + *¥*;*j) 


(B.31) 


Equation (B.28) can be written as 




= 0 


(B.32) 


To facilitate the application of approximate factorization, the cross-derivative terms 
are lagged in time in the usual way to obtain the relaxation algorithm 


.A'*' 


,A^ . 


1 / - a ( ](^ +1 - r ) 

II W 1 _ II 

/. w {« 5 (^r + 2( T )" + , + iA ? |i ? ( 7 ri} 


(B.33) 


Equation (B.32) is of the same form as equation (B.16) and hence tangency can 
be enforced implicitly in the ADI algorithm given by equation (B.19) with the 
tridiagonal and right-hand-side terms modified appropriately. 
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B. Boundary Layer Equations 

U Governing equations 

In two dimensions, the unsteady, incompressible laminar boundary-layer equa- 
tions are 

u t -r - 0 (B.34) 

pu t -f p{uu.j. 4- iviiz) - ~px Uzz (B.35) 

In these equations, the normal coordinate z and the normal component of velocity 
w are scaled by the square root of the free stream Reynolds number Re^ . 

Using a general coordinate transformation defined by 


c — 

S — 


C(r\. 




A — 


( 

v~ 


•y\ 

~r 


the above equations are transformed as 


Uc£x + Uc$x + w^z = 0 


pu t + p{Uu£ + W Uf) = ~p^x + (Wc&K-Ci: 


where U and W are unsealed contravariant velocities defined as 


U — & + £x u 


w = ft + + fete 


' ' / 


(B.37) 

(B.38) 


(B.39) 


2, Numerical algorithms 

The boundary-layer equations are solved with a relaxation algorithm which can 
be run in either a time-accurate or an iterative mode. The relaxation algorithm 
was designed by Van Dalsem and Steger (1985) to yield convergence to a steady- 
state algorithm quickly. The continuity and momentum equations are solved in an 
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uncoupled manner at each time step. The momentum equation is used to update 
the streamwise velocity u using upwind differencing in the streamwise direction and 
central differencing in the normal direction. The continuity equation is used to 
update the normal velocity w using second-order-accurate central differencing in 
the ^-direction. The trapezoidal rule is used to integrate the continuity equation in 
the f-direction to obtain w. 

The algorithm was designed for three-dimensional unsteady compressible flows 
by Van Dalsem and Steger (1985). A simplified version for two-dimensional incom- 
pressible flows which was used in this study is presented below. 

1) Update u at the new time step from the x-momentum equation 

pV<ii n+1 ■+■ p(U + 6cu n * 1 + U~ d{u n + IV 6 s *u n+1 ) = (B.40) 

2) Integrate the continuity equation for w using updated values of ii 

V ? 1 B n+1 = UjScil + fe$ ? u) n+I (B.41) 

3) Repeat steps (l)-(2) till convergence is obtained to the steady-state solution. 



